TABLE OF CONTENTS
- ABINIT/greendftcompute_green
- ABINIT/m_green
- m_green/add_int_fct
- m_green/check_fourier_green
- m_green/compa_occup_ks
- m_green/compute_green
- m_green/copy_green
- m_green/destroy_green
- m_green/destroy_green_tau
- m_green/distrib_paral
- m_green/fermi_green
- m_green/fourier_fct
- m_green/fourier_green
- m_green/green_type
- m_green/icip_green
- m_green/init_green
- m_green/init_green_tau
- m_green/int_fct
- m_green/integrate_green
- m_green/local_ks_green
- m_green/newton
- m_green/occup_green_tau
- m_green/occupfd
- m_green/print_green
- m_green/printocc_green
- m_green/spline_fct
- newton/compute_nb_elec
- newton/function_and_deriv
ABINIT/greendftcompute_green [ Functions ]
NAME
greendftcompute_green
FUNCTION
Compute levels for ctqmc
COPYRIGHT
Copyright (C) 1999-2024 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
SOURCE
2895 subroutine greendftcompute_green(cryst_struc,green,pawang,paw_dmft) 2896 2897 !Arguments ------------------------------------ 2898 !scalars 2899 type(crystal_t),intent(in) :: cryst_struc 2900 type(pawang_type), intent(in) :: pawang 2901 type(paw_dmft_type), intent(in) :: paw_dmft 2902 type(green_type),intent(inout) :: green 2903 2904 !Local variables ------------------------------ 2905 ! scalars 2906 integer :: ifreq,iband,ikpt,isppol 2907 integer :: mbandc,natom,nspinor,nsppol,nkpt 2908 character(len=500) :: message 2909 ! arrays 2910 !************************************************************************ 2911 2912 mbandc=paw_dmft%mbandc 2913 nkpt=paw_dmft%nkpt 2914 nsppol=paw_dmft%nsppol 2915 nspinor=paw_dmft%nspinor 2916 natom=paw_dmft%natom 2917 2918 ! write(6,*) green%oper(1)%ks(1,1,1,1) 2919 if(green%oper(1)%has_operks==0) then 2920 ABI_ERROR("greendft%oper(1)%ks not allocated") 2921 endif 2922 2923 !====================================== 2924 !Get Green's function G=1/(iw_n+mu-e_nk) 2925 !====================================== 2926 do ifreq=1,green%nw 2927 do iband=1,mbandc 2928 do ikpt=1,nkpt 2929 do isppol=1,nsppol 2930 ! write(6,*) green%oper(ifreq)%ks(isppol,ikpt,iband,iband) 2931 ! write(6,*) ifreq,iband,ikpt,isppol 2932 ! write(6,*) cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp) 2933 ! write(6,*) paw_dmft%fermie 2934 ! write(6,*) paw_dmft%eigen_dft(isppol,ikpt,iband) 2935 green%oper(ifreq)%ks(isppol,ikpt,iband,iband)=& 2936 cone/(cmplx(0.d0,paw_dmft%omega_lo(ifreq),kind=dp)+paw_dmft%fermie-paw_dmft%eigen_dft(isppol,ikpt,iband)) 2937 end do 2938 end do 2939 end do 2940 2941 2942 !====================================================================== 2943 ! Compute local Green's function 2944 !====================================================================== 2945 call loc_oper(green%oper(ifreq),paw_dmft,1) 2946 2947 !====================================================================== 2948 ! Symetrize 2949 !====================================================================== 2950 call sym_matlu(cryst_struc,green%oper(ifreq)%matlu,pawang,paw_dmft) 2951 enddo 2952 write(message,'(a,2x,a,f13.5)') ch10," == Print DFT Green's function for last frequency" 2953 call wrtout(std_out,message,'COLL') 2954 call print_matlu(green%oper(paw_dmft%dmft_nwlo)%matlu,natom,1) 2955 2956 end subroutine greendftcompute_green
ABINIT/m_green [ Modules ]
NAME
m_green
FUNCTION
COPYRIGHT
Copyright (C) 2006-2024 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 .
SOURCE
15 #if defined HAVE_CONFIG_H 16 #include "config.h" 17 #endif 18 19 20 #include "abi_common.h" 21 22 MODULE m_green 23 24 use defs_basis 25 use m_xmpi 26 use m_abicore 27 use m_errors 28 use m_lib_four 29 use m_splines 30 31 !use defs_abitypes, only : MPI_type 32 use m_crystal, only : crystal_t 33 use m_oper, only : init_oper, destroy_oper, copy_oper, print_oper, inverse_oper, upfold_oper, & 34 loc_oper, trace_oper, oper_type 35 use m_paw_dmft, only: paw_dmft_type, construct_nwli_dmft 36 use m_matlu, only : diff_matlu,print_matlu,trace_matlu, matlu_type, & 37 sym_matlu,zero_matlu,add_matlu,shift_matlu,init_matlu,destroy_matlu,copy_matlu 38 use m_fstrings, only : int2char4 39 use m_pawang, only : pawang_type 40 use m_self, only : self_type 41 use m_io_tools, only : open_file 42 use m_time, only : timab 43 44 implicit none 45 46 private 47 48 public :: init_green 49 public :: destroy_green 50 public :: init_green_tau 51 public :: destroy_green_tau 52 public :: print_green 53 public :: printocc_green 54 public :: compute_green 55 public :: integrate_green 56 public :: icip_green 57 public :: fourier_green 58 public :: check_fourier_green 59 public :: compa_occup_ks 60 public :: copy_green 61 public :: occup_green_tau 62 public :: add_int_fct 63 public :: int_fct 64 public :: fourier_fct 65 public :: spline_fct 66 private :: occupfd 67 private :: distrib_paral 68 public :: greendftcompute_green 69 public :: fermi_green 70 public :: newton 71 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-2024 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
SOURCE
2296 subroutine add_int_fct(ifreq,ff,ldiag,omega_current,option,integral,temp,wgt_wlo,dmft_nwlo) 2297 2298 !Arguments ------------------------------------ 2299 !type 2300 integer,intent(in) :: ifreq 2301 logical,intent(in) :: ldiag 2302 integer,intent(in) :: option,dmft_nwlo 2303 complex(dpc),intent(inout) :: integral 2304 complex(dpc), intent(in) :: ff 2305 complex(dpc), intent(in) :: omega_current 2306 real(dp), intent(in) :: temp, wgt_wlo 2307 2308 !local variables------------------------------- 2309 real(dp) :: omega 2310 ! ********************************************************************* 2311 omega=aimag(omega_current) 2312 if(ldiag) then 2313 2314 if(option==1) then ! nspinor==1 2315 ! write(500,*) paw_dmft%omega_lo(ifreq),real(ff(ifreq)),imag(ff(ifreq)) 2316 integral=integral+2.d0*temp * & 2317 & real( ff-one / ( j_dpc*omega ) ) * & 2318 & wgt_wlo 2319 if(ifreq==dmft_nwlo) integral=integral+half 2320 ! integral=integral+half 2321 ! the if is here, to count only one time this correction 2322 endif 2323 2324 if(option==2) then ! nspinor==2 2325 integral=integral+2.d0*temp * & 2326 & ( ff-one / ( j_dpc*omega ) ) * & 2327 & wgt_wlo 2328 if(ifreq==dmft_nwlo) integral=integral+half 2329 ! integral=integral+half 2330 ! the if is here, to count only one time this correction 2331 endif 2332 2333 2334 else ! ldiag 2335 2336 ! write(std_out,*) "nondiag" 2337 if(option==1) then 2338 integral=integral+2.d0*temp * & 2339 & real( ff ) * & 2340 & wgt_wlo 2341 endif 2342 if(option==2) then 2343 integral=integral+2.d0*temp * & 2344 & ff * & 2345 & wgt_wlo 2346 endif 2347 endif ! ldiag 2348 2349 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
SOURCE
2148 subroutine check_fourier_green(cryst_struc,green,paw_dmft,pawang) 2149 2150 !Arguments ------------------------------------ 2151 !type 2152 type(crystal_t),intent(in) :: cryst_struc 2153 type(green_type),intent(inout) :: green 2154 !type(MPI_type), intent(in) :: mpi_enreg 2155 type(pawang_type),intent(in) :: pawang 2156 type(paw_dmft_type), intent(inout) :: paw_dmft 2157 2158 !local variables------------------------------- 2159 type(green_type) :: green_check 2160 character(len=500) :: message 2161 ! ********************************************************************* 2162 ! Only imaginary frequencies here 2163 if(green%w_type=="real") then 2164 message = 'check_fourier_green not implemented for real frequency' 2165 ABI_BUG(message) 2166 endif 2167 2168 call init_green(green_check,paw_dmft) 2169 call init_green_tau(green_check,paw_dmft) 2170 call copy_green(green,green_check,opt_tw=2) 2171 2172 write(message,'(2a,i3,13x,a)') ch10,' === Inverse Fourier Transform w->t of Weiss Field' 2173 call wrtout(std_out,message,'COLL') 2174 call fourier_green(cryst_struc,green_check,paw_dmft,& 2175 & pawang,opt_ksloc=2,opt_tw=-1) 2176 2177 write(message,'(3a)') ch10,' === Print (for check by user) of occupation matrix'& 2178 & ,' after fourier transform with respect to initial one' 2179 call wrtout(std_out,message,'COLL') 2180 call printocc_green(green_check,6,paw_dmft,3) 2181 2182 write(message,'(2a,i3,13x,a)') ch10,' === Direct Fourier Transform t->w of Weiss Field' 2183 call wrtout(std_out,message,'COLL') 2184 call fourier_green(cryst_struc,green_check,paw_dmft,& 2185 & pawang,opt_ksloc=2,opt_tw=1) 2186 ! call print_matlu(green%oper(1)%matlu,paw_dmft%natom,1) ! debug 2187 2188 call integrate_green(cryst_struc,green_check,paw_dmft,& 2189 & pawang,prtopt=2,opt_ksloc=2) 2190 2191 write(message,'(3a)') ch10,' === Print (for check by user) of occupation matrix'& 2192 & ,' after double fourier transform with respect to initial one' 2193 call wrtout(std_out,message,'COLL') 2194 call printocc_green(green_check,5,paw_dmft,3) 2195 2196 call destroy_green_tau(green_check) 2197 call destroy_green(green_check) 2198 2199 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
SOURCE
2217 subroutine compa_occup_ks(green,paw_dmft) 2218 2219 !Arguments ------------------------------------ 2220 !type 2221 type(green_type),intent(inout) :: green 2222 type(paw_dmft_type), intent(inout) :: paw_dmft 2223 2224 !local variables------------------------------- 2225 character(len=500) :: message 2226 integer :: ib,ikpt,isppol 2227 integer :: ib_,ikpt_,isppol_ 2228 integer :: ib__,ikpt__,isppol__ 2229 real(dp) :: diffmax,occ1,occ2,occ_a,occ_b,diffrel,occ_aa,occ_bb 2230 ! ********************************************************************* 2231 diffmax=zero 2232 diffrel=zero 2233 do isppol=1,paw_dmft%nsppol 2234 do ikpt=1,paw_dmft%nkpt 2235 do ib=1,paw_dmft%mbandc 2236 occ1=occupfd(paw_dmft%eigen_dft(isppol,ikpt,ib),paw_dmft%fermie,paw_dmft%temp) 2237 occ2=real(green%occup%ks(isppol,ikpt,ib,ib)) 2238 if(abs(occ1-occ2)>diffmax) then 2239 diffmax=abs(occ1-occ2) 2240 occ_a=occ1 2241 occ_b=occ2 2242 ib_=ib;isppol_=isppol;ikpt_=ikpt 2243 endif 2244 if(abs(two*(occ1-occ2)/(occ1+occ2))>diffrel) then 2245 diffrel=abs(two*(occ1-occ2)/(occ1+occ2)) 2246 occ_aa=occ1 2247 occ_bb=occ2 2248 ib__=ib;isppol__=isppol;ikpt__=ikpt 2249 endif 2250 enddo 2251 enddo 2252 enddo 2253 2254 2255 write(message,'(2a)') ch10,' === Compare green function occupations and Fermi Dirac occupations' 2256 call wrtout(std_out,message,'COLL') 2257 write(message,'(2a,f12.5,2a,f12.5,2a,f12.5,2a,3i5)') ch10,' = Max difference is',diffmax,& 2258 & ch10,' Corresponding Occupation from green function is',occ_b,& 2259 & ch10,' Corresponding Occupation Fermi Dirac weight is',occ_a,& 2260 & ch10,' (For polarization, k-point and band index) ',isppol_,ikpt_,ib_ 2261 call wrtout(std_out,message,'COLL') 2262 write(message,'(2a,f12.5,2a,f12.5,2a,f12.5,2a,3i5)') ch10,' = Max relative difference is',diffrel,& 2263 & ch10,' Corresponding Occupation from green function is',occ_bb,& 2264 & ch10,' Corresponding Occupation Fermi Dirac weight is',occ_aa,& 2265 & ch10,' (For polarization, k-point and band index) ',isppol__,ikpt__,ib__ 2266 call wrtout(std_out,message,'COLL') 2267 2268 end subroutine compa_occup_ks
m_green/compute_green [ Functions ]
[ Top ] [ m_green ] [ Functions ]
NAME
compute_green
FUNCTION
compute green function from DFT 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
SOURCE
956 subroutine compute_green(cryst_struc,green,paw_dmft,pawang,prtopt,self,opt_self,& 957 & opt_nonxsum,opt_nonxsum2) 958 959 !Arguments ------------------------------------ 960 !type 961 type(crystal_t),intent(in) :: cryst_struc 962 type(green_type),intent(inout) :: green 963 type(paw_dmft_type), intent(in) :: paw_dmft 964 !type(MPI_type), intent(in) :: mpi_enreg 965 type(pawang_type),intent(in) :: pawang 966 type(self_type), intent(inout) :: self 967 integer, intent(in) :: prtopt 968 integer, optional, intent(in) :: opt_self 969 integer, optional, intent(in) :: opt_nonxsum 970 integer, optional, intent(in) :: opt_nonxsum2 971 972 !local variables------------------------------- 973 !logical :: lintegrate 974 integer :: iatom,ib,ib1,ierr,ifreq,ikpt,is 975 integer :: icomp_chloc 976 integer :: natom,mband,mbandc,nkpt 977 integer :: myproc,nproc,spacecomm,nspinor,nsppol,option,optself,optnonxsum 978 integer :: optnonxsum2 979 integer, allocatable :: procb_ifreq(:) 980 character(len=500) :: message 981 real(dp) :: fermilevel 982 real(dp), allocatable :: Id(:,:) 983 ! integer, allocatable :: procb(:,:),proct(:,:) 984 real(dp) :: tsec(2) 985 type(oper_type) :: self_minus_hdc_oper 986 complex(dpc) :: omega_current 987 type(oper_type) :: green_temp 988 ! ********************************************************************* 989 990 !lintegrate=.true. 991 !if(lintegrate.and.green%w_type=="real") then 992 !if(green%w_type=="real") then 993 ! message = 'integrate_green not implemented for real frequency' 994 ! ABI_BUG(message) 995 !endif 996 call timab(624,1,tsec) 997 if(present(opt_self)) then 998 optself=opt_self 999 else 1000 optself=0 1001 endif 1002 if(present(opt_nonxsum)) then 1003 optnonxsum=opt_nonxsum 1004 else 1005 optnonxsum=0 1006 endif 1007 if(present(opt_nonxsum2)) then 1008 optnonxsum2=opt_nonxsum2 1009 else 1010 optnonxsum2=0 1011 endif 1012 1013 if(prtopt>0) then 1014 write(message,'(2a,i3,13x,a)') ch10,' === Compute green function ' 1015 call wrtout(std_out,message,'COLL') 1016 endif 1017 1018 if(self%nw/=green%nw) then 1019 message = ' BUG: frequencies for green and self not coherent' 1020 ABI_BUG(message) 1021 endif 1022 1023 ! Initialise spaceComm, myproc, and nproc 1024 spacecomm=paw_dmft%spacecomm 1025 myproc=paw_dmft%myproc 1026 nproc=paw_dmft%nproc 1027 1028 ! Initialise integers 1029 mband = paw_dmft%mband 1030 mbandc = paw_dmft%mbandc 1031 nkpt = paw_dmft%nkpt 1032 nspinor = paw_dmft%nspinor 1033 nsppol = paw_dmft%nsppol 1034 natom = paw_dmft%natom 1035 1036 ! Initialise for compiler 1037 omega_current=czero 1038 icomp_chloc=0 1039 1040 ! ============================== 1041 ! Initialise Identity 1042 ! ============================== 1043 ABI_MALLOC(Id,(paw_dmft%mbandc,paw_dmft%mbandc)) 1044 Id=zero 1045 do ib=1, paw_dmft%mbandc 1046 Id(ib,ib)=one 1047 enddo ! ib 1048 1049 if(prtopt/=0.and.prtopt>-100) then 1050 write(message,'(2a)') ch10,' == Green function is computed:' 1051 call wrtout(std_out,message,'COLL') 1052 endif 1053 option=1 1054 fermilevel=paw_dmft%fermie 1055 if(option==123) then 1056 fermilevel=2.d0 1057 write(message,'(2a,e14.3,a)') ch10,& 1058 & ' Warning (special case for check: fermi level=',fermilevel,')' 1059 call wrtout(std_out,message,'COLL') 1060 endif 1061 1062 ! ==================================================== 1063 ! Upfold self-energy and double counting Self_imp -> self(k) 1064 ! ==================================================== 1065 ! if(optself==1) then 1066 ! do ifreq=1,green%nw 1067 ! call upfold_oper(self%oper(ifreq),paw_dmft,1) 1068 ! enddo ! ifreq 1069 ! call upfold_oper(self%hdc,paw_dmft,1) 1070 ! endif 1071 1072 ! ================================================================= 1073 ! Initialize green%oper before calculation (important for xmpi_sum) 1074 ! ================================================================= 1075 do ifreq=1,green%nw 1076 call zero_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom) 1077 enddo ! ifreq 1078 !call xmpi_barrier(spacecomm) 1079 1080 ! ================================ 1081 ! == Compute Green function G(k) 1082 ! ================================ 1083 green%occup%ks=czero 1084 ABI_MALLOC(procb_ifreq,(paw_dmft%nkpt)) 1085 do ifreq=1,green%nw 1086 !if(present(iii)) write(6,*) ch10,'ifreq self', ifreq,self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) 1087 ! ==================================================== 1088 ! First Upfold self-energy and double counting Self_imp -> self(k) 1089 ! ==================================================== 1090 ! if(mod(ifreq-1,nproc)==myproc) then 1091 ! write(6,*) "compute_green ifreq",ifreq, mod(ifreq-1,nproc)==myproc,proct(ifreq,myproc)==1 1092 if(green%proct(ifreq,myproc)==1) then 1093 if(green%w_type=="imag") then 1094 omega_current=cmplx(zero,green%omega(ifreq),kind=dp) 1095 else if(green%w_type=="real") then 1096 omega_current=cmplx(green%omega(ifreq),paw_dmft%temp,kind=dp) 1097 endif 1098 call init_oper(paw_dmft,self_minus_hdc_oper) 1099 call init_oper(paw_dmft,green_temp) 1100 1101 call add_matlu(self%oper(ifreq)%matlu,self%hdc%matlu,& 1102 & self_minus_hdc_oper%matlu,green%oper(ifreq)%natom,-1) 1103 ! do iatom = 1 , natom 1104 !write(6,*) 'self matlu', ifreq, self%oper(ifreq)%matlu(1)%mat(1,1,1,1,1) 1105 !write(6,*) 'self hdc ', ifreq, self%hdc%matlu(1)%mat(1,1,1,1,1) 1106 !write(6,*) 'self_minus_hdc_oper ', ifreq, self_minus_hdc_oper%matlu(1)%mat(1,1,1,1,1) 1107 ! enddo ! natom 1108 if(paw_dmft%dmft_solv==4) then 1109 call shift_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom,cmplx(self%qmc_shift,0.d0,kind=dp),-1) 1110 call shift_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom,cmplx(self%qmc_xmu,0.d0,kind=dp),-1) 1111 endif 1112 if(optself==0) then 1113 call zero_matlu(self_minus_hdc_oper%matlu,paw_dmft%natom) 1114 end if 1115 1116 procb_ifreq=green%procb(ifreq,:) 1117 call upfold_oper(self_minus_hdc_oper,paw_dmft,1,procb=procb_ifreq,iproc=myproc,prt=1) 1118 do ib1 = 1 , paw_dmft%mbandc 1119 do ib = 1 , paw_dmft%mbandc 1120 do ikpt = 1 , paw_dmft%nkpt 1121 do is = 1 , paw_dmft%nsppol 1122 if (green%procb(ifreq,ikpt)==myproc) then 1123 ! green%oper(ifreq)%ks(is,ikpt,ib,ib1)= & 1124 green_temp%ks(is,ikpt,ib,ib1)= & 1125 & ( omega_current & 1126 & + fermilevel & 1127 & - paw_dmft%eigen_dft(is,ikpt,ib)) * Id(ib,ib1) & 1128 & - self_minus_hdc_oper%ks(is,ikpt,ib,ib1) 1129 !if(ikpt==2.and.ib==ib1) then 1130 ! write(6,*) "self",ib1,ib,ikpt,is,ifreq,self_minus_hdc_oper%ks(is,ikpt,ib,ib1) 1131 !endif 1132 !& - (self%oper(ifreq)%ks(is,ikpt,ib,ib1)-self%hdc%ks(is,ikpt,ib,ib1)) 1133 ! if(prtopt>5) then 1134 ! if(ikpt==1.and.(ifreq==1.or.ifreq==3).and.ib==16.and.ib1==16) then 1135 ! write(std_out,*) 'omega_current ',omega_current 1136 ! write(std_out,*) 'fermilevel ',fermilevel 1137 ! write(std_out,*) ' paw_dmft%eigen_dft(is,ikpt,ib) ', paw_dmft%eigen_dft(is,ikpt,ib),Id(ib,ib1) 1138 ! write(std_out,*) 'self_minus_hdc_oper%ks(is,ikpt,ib,ib1)',self_minus_hdc_oper%ks(is,ikpt,ib,ib1) 1139 ! write(std_out,*) 'green ',green%oper(ifreq)%ks(is,ikpt,ib,ib1) 1140 ! endif 1141 ! if(ib==1.and.ib1==3) then 1142 ! write(std_out,*) "ff compute",ikpt,ifreq,is,ikpt,ib,ib1 1143 ! write(std_out,*) "ff compute",ikpt,ifreq, green_temp%ks(is,ikpt,ib,ib1) 1144 ! write(std_out,*) "ff details",paw_dmft%eigen_dft(is,ikpt,ib) 1145 ! write(std_out,*) "ff details2",fermilevel 1146 ! write(std_out,*) "ff details3",Id(ib,ib1) 1147 ! ! write(std_out,*) "ff details4",self_minus_hdc_oper%ks(is,ikpt,ib,ib1) 1148 ! endif 1149 endif 1150 enddo ! is 1151 enddo ! ikpt 1152 enddo ! ib 1153 enddo ! ib1 1154 ! call print_oper(green%oper(ifreq),9,paw_dmft,3) 1155 ! write(std_out,*) 'after print_oper' 1156 ! if(ifreq==1.or.ifreq==3) then 1157 ! write(std_out,*) 'after print_oper', ifreq 1158 ! write(std_out,*) 'green1 ifreq %ks(1,1,16,16)',ifreq,green%oper(ifreq)%ks(1,1,16,16) 1159 ! endif 1160 ! write(std_out,*) 'before inverse_oper' 1161 call inverse_oper(green_temp,2,prtopt,procb=procb_ifreq,iproc=myproc) 1162 if (green%oper(1)%has_operks==1) green%oper(ifreq)%ks=green_temp%ks 1163 !if(ifreq==1) then 1164 ! write(std_out,*) "1188",green_temp%ks(1,1,8,8) 1165 ! write(std_out,*) "1189",green_temp%ks(1,1,8,9) 1166 ! write(std_out,*) "1198",green_temp%ks(1,1,9,8) 1167 ! write(std_out,*) "1199",green_temp%ks(1,1,9,9) 1168 !endif 1169 1170 !if(lintegrate) then 1171 ! accumulate integration 1172 if(green%w_type/="real") then 1173 do is = 1 , paw_dmft%nsppol 1174 do ib = 1 , paw_dmft%mbandc 1175 do ib1 = 1 , paw_dmft%mbandc 1176 do ikpt = 1 , paw_dmft%nkpt 1177 if (green%procb(ifreq,ikpt)==myproc) then 1178 call add_int_fct(ifreq,green_temp%ks(is,ikpt,ib,ib1),ib==ib1, & 1179 & omega_current,2,green%occup%ks(is,ikpt,ib,ib1), & 1180 & paw_dmft%temp,paw_dmft%wgt_wlo(ifreq),paw_dmft%dmft_nwlo) 1181 endif 1182 enddo ! ikpt 1183 enddo ! ib1 1184 enddo ! ib 1185 enddo ! is 1186 endif 1187 ! write(std_out,*) 'after inverse_oper' 1188 ! if(ikpt==1.and.is==1.and.ib==1.and.ib1==1) then 1189 ! write(6,*) 'occup(is,ikpt,ib,ib1)',ifreq,green%occup%ks(1,1,1,1),green_temp%ks(1,1,1,1) 1190 ! endif 1191 ! write(std_out,*) 'green1afterinversion %ks(1,1,16,16)',ifreq,green%oper(ifreq)%ks(1,1,16,16) 1192 ! endif 1193 ! write(std_out,*) 'before flush' 1194 ! call flush(std_out) 1195 !endif 1196 ! ================================ 1197 ! == Compute Local Green function 1198 ! ================================ 1199 !write(message,'(2a)') ch10,' loc' 1200 !call wrtout(std_out,message,'COLL') 1201 !call flush(std_out) 1202 call loc_oper(green_temp,paw_dmft,1,procb=procb_ifreq,iproc=myproc) 1203 !if(ifreq==1) then 1204 ! write(std_out,*) "4411",green_temp%matlu(1)%mat(4,4,1,1,1) 1205 ! write(std_out,*) "4512",green_temp%matlu(1)%mat(4,5,1,1,2) 1206 ! write(std_out,*) "5421",green_temp%matlu(1)%mat(5,4,1,2,1) 1207 ! write(std_out,*) "5522",green_temp%matlu(1)%mat(5,5,1,2,2) 1208 ! write(std_out,*) "(5512)",green_temp%matlu(1)%mat(5,5,1,1,2) 1209 !endif 1210 ! write(std_out,*) ifreq,nproc,'before if after loc_oper' 1211 ! if(ifreq==1.or.ifreq==11) then 1212 ! write(std_out,*) ifreq,nproc,'before sym' 1213 ! call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0) 1214 ! ! ok 1215 ! endif 1216 ! call flush(std_out) 1217 call sym_matlu(cryst_struc,green_temp%matlu,pawang,paw_dmft) 1218 call copy_matlu(green_temp%matlu,green%oper(ifreq)%matlu,natom) 1219 ! if(ifreq==1.and.ifreq==11) then 1220 ! write(std_out,*) ifreq,nproc,'after sym' 1221 ! call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0) 1222 ! ! ok 1223 ! endif 1224 call destroy_oper(self_minus_hdc_oper) 1225 call destroy_oper(green_temp) 1226 ! call flush(std_out) 1227 endif ! parallelisation 1228 enddo ! ifreq 1229 ABI_FREE(procb_ifreq) 1230 1231 ! ============================================= 1232 ! built total green function (sum over procs). 1233 ! ============================================= 1234 !call xmpi_barrier(spacecomm) 1235 do ifreq=1,green%nw 1236 ! call xmpi_sum(green%oper(ifreq)%ks,spacecomm,ierr) 1237 ! print *, "myproc, proct, ifreq ------------------- ", myproc, ifreq 1238 ! do ikpt=1,paw_dmft%nkpt 1239 ! call xmpi_bcast(green%oper(ifreq)%ks(:,ikpt,:,:),procb(ifreq,ikpt),spacecomm,ierr) 1240 ! enddo 1241 !! or 1242 if(optnonxsum==0.and.green%oper(ifreq)%has_operks==1) then 1243 call xmpi_sum(green%oper(ifreq)%ks,spacecomm,ierr) 1244 else if(optnonxsum==0.and.green%oper(ifreq)%has_operks==0) then 1245 message = 'optnonxsum==0 and green%oper(ifreq)%has_operks==0: not compatible' 1246 ABI_BUG(message) 1247 endif 1248 ! endif 1249 ! enddo ! ifreq 1250 ! print *,"myproc", myproc 1251 if(optnonxsum2==0) then 1252 do iatom=1,green%oper(ifreq)%natom 1253 if(green%oper(ifreq)%matlu(iatom)%lpawu.ne.-1) then 1254 call xmpi_sum(green%oper(ifreq)%matlu(iatom)%mat,spacecomm,ierr) 1255 endif 1256 enddo ! iatom 1257 green%has_greenmatlu_xsum=1 1258 else if(optnonxsum2==1) then 1259 green%has_greenmatlu_xsum=0 1260 endif 1261 ! if(ifreq==1.or.ifreq==11) then 1262 ! write(std_out,*) ifreq,nproc,'after xsum' 1263 ! call print_matlu(green%oper(ifreq)%matlu,green%oper(ifreq)%natom,2,-1,0) 1264 ! ! ok 1265 ! endif 1266 enddo ! ifreq 1267 !call xmpi_barrier(spacecomm) 1268 call xmpi_sum(green%occup%ks,spacecomm,ierr) 1269 ! 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) 1270 1271 if(prtopt/=0.and.prtopt>-100) then 1272 write(message,'(2a)') ch10,& 1273 & ' == Local Green function has been computed and projected on local orbitals' 1274 call wrtout(std_out,message,'COLL') 1275 endif 1276 ! useless test 1277 if(abs(prtopt)>=4.and.prtopt>-100) then 1278 write(message,'(2a)') ch10,' == Green function is now printed for first frequency' 1279 call wrtout(std_out,message,'COLL') 1280 call print_oper(green%oper(1),9,paw_dmft,3) 1281 write(message,'(2a)') ch10,' == Green function is now printed for second frequency' 1282 call wrtout(std_out,message,'COLL') 1283 call print_oper(green%oper(2),9,paw_dmft,3) 1284 if(paw_dmft%dmft_nwlo>=11) then 1285 write(message,'(2a)') ch10,' == Green function is now printed for 11th frequency' 1286 call wrtout(std_out,message,'COLL') 1287 call print_oper(green%oper(11),9,paw_dmft,3) 1288 endif 1289 endif 1290 ! call flush(std_out) 1291 1292 ABI_FREE(Id) 1293 call timab(624,2,tsec) 1294 1295 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
SOURCE
479 subroutine copy_green(green1,green2,opt_tw) 480 481 !Arguments ------------------------------------ 482 !type 483 type(green_type),intent(in) :: green1 484 type(green_type),intent(inout) :: green2 485 integer, intent(in) :: opt_tw 486 487 !local variables------------------------------- 488 integer :: ifreq,itau 489 ! ********************************************************************* 490 491 if(opt_tw==2) then 492 call copy_oper(green1%occup,green2%occup) 493 do ifreq=1, green1%nw 494 call copy_oper(green1%oper(ifreq),green2%oper(ifreq)) 495 if(green1%has_greenmatlu_xsum==1) then ! indicate to integrate_green that xsum has been done 496 ! for matlu in compute_green. 497 green2%has_greenmatlu_xsum=1 498 endif 499 enddo 500 else if (opt_tw==1) then 501 call copy_oper(green1%occup_tau,green2%occup_tau) 502 do itau=1,green1%dmftqmc_l 503 call copy_oper(green1%oper_tau(itau),green2%oper_tau(itau)) 504 enddo 505 endif 506 507 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
SOURCE
371 subroutine destroy_green(green) 372 373 !Arguments ------------------------------------ 374 !scalars 375 type(green_type),intent(inout) :: green 376 377 !local variables------------------------------- 378 integer :: ifreq 379 380 ! ********************************************************************* 381 call destroy_oper(green%occup) 382 if ( allocated(green%oper)) then 383 do ifreq=1,green%nw 384 call destroy_oper(green%oper(ifreq)) 385 enddo 386 ABI_FREE(green%oper) 387 end if 388 if ( allocated(green%charge_matlu)) then 389 ABI_FREE(green%charge_matlu) 390 end if 391 green%has_charge_matlu=0 392 393 if ( allocated(green%charge_matlu_prev)) then 394 ABI_FREE(green%charge_matlu_prev) 395 end if 396 green%has_charge_matlu_prev=0 397 398 if ( allocated(green%charge_matlu_solver)) then 399 ABI_FREE(green%charge_matlu_solver) 400 end if 401 green%has_charge_matlu_solver=0 402 403 if ( allocated(green%ecorr_qmc)) then 404 ABI_FREE(green%ecorr_qmc) 405 end if 406 if ( allocated(green%procb)) then 407 ABI_FREE(green%procb) 408 end if 409 if ( allocated(green%proct)) then 410 ABI_FREE(green%proct) 411 end if 412 green%omega => null() 413 414 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
SOURCE
431 subroutine destroy_green_tau(green) 432 433 !Arguments ------------------------------------ 434 !scalars 435 type(green_type),intent(inout) :: green 436 ! integer, optional, intent(in) :: opt_ksloc 437 !local variables------------------------------- 438 integer :: itau 439 ! integer :: optksloc 440 ! ********************************************************************* 441 ! if(present(opt_ksloc)) then 442 ! optksloc=opt_ksloc 443 ! else 444 ! optksloc=3 445 ! endif 446 447 call destroy_oper(green%occup_tau) 448 if ( allocated(green%oper_tau)) then 449 do itau=1,green%dmftqmc_l 450 call destroy_oper(green%oper_tau(itau)) 451 enddo 452 ABI_FREE(green%oper_tau) 453 end if 454 if ( allocated(green%tau)) then 455 ABI_FREE(green%tau) 456 end if 457 458 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.
SOURCE
2777 subroutine distrib_paral(nkpt,nproc,nw,nw_perproc,procb,proct) 2778 2779 !Arguments ------------------------------------ 2780 !type 2781 ! Integrate analytic tail 1/(iw-mu) 2782 integer, intent(in) :: nw,nkpt,nproc 2783 integer, intent(out) :: nw_perproc 2784 integer, intent(out):: procb(nw,nkpt),proct(nw,0:nproc-1) 2785 !Local variables------------------------------- 2786 integer :: ikpt,iw,ir,ratio,m,iproc 2787 integer, allocatable :: proca(:,:) 2788 ! ********************************************************************* 2789 2790 2791 proct(:,:)=0 2792 procb(:,:)=-1 2793 2794 if(nproc.ge.2*nw) then 2795 !write(6,*) "AA" 2796 ratio=nproc/nw 2797 ABI_MALLOC(proca,(nw,nproc/nw)) 2798 do ir=1,ratio 2799 do iw=1,nw 2800 proca(iw,ir)=iw+(ir-1)*nw-1 2801 proct(iw,iw+(ir-1)*nw-1)=1 2802 enddo 2803 enddo 2804 ! do iw=1,nw 2805 ! write(6,*) " freq, procs" 2806 ! write(6,*) iw,proca(iw,:) 2807 ! enddo 2808 do iw=1,nw 2809 do ikpt=1,nkpt 2810 if(ikpt.le.ratio) procb(iw,ikpt)=proca(iw,ikpt) 2811 if(ikpt.ge.ratio) then 2812 m=mod(ikpt,ratio) 2813 if(m==0) then 2814 procb(iw,ikpt)=proca(iw,ratio) 2815 else 2816 procb(iw,ikpt)=proca(iw,m) 2817 endif 2818 endif 2819 ! write(6,*) " freq, k-point, proc" 2820 ! write(6,*) iw,ikpt,procb(iw,ikpt) 2821 enddo 2822 enddo 2823 nw_perproc=1 2824 ABI_FREE(proca) 2825 2826 else if (nproc.ge.nw) then 2827 !write(6,*) "BB" 2828 do iw=1,nw 2829 procb(iw,:)= iw 2830 proct(iw,iw)=1 2831 enddo 2832 ! do iw=1,nw 2833 ! write(6,*) " freq, proc" 2834 ! write(6,*) iw, procb(iw,1) 2835 ! enddo 2836 nw_perproc=1 2837 2838 else if (nproc.le.nw) then 2839 !write(6,*) "CC" 2840 ratio=nw/nproc 2841 !write(6,*) "ratio", ratio 2842 do iproc=0,nproc-1 2843 do iw=1,nw 2844 if (mod(iw-1,nproc)==iproc) then 2845 procb(iw,:)=iproc 2846 proct(iw,iproc)=1 2847 endif 2848 enddo 2849 enddo 2850 ! do iw=1,nw 2851 ! write(6,*) "iw, iproc", iw, procb(iw,1) 2852 ! enddo 2853 nw_perproc=ratio+1 2854 ! some procs will compute a number of frequency which is ratio and some 2855 ! other will compute a number of frequency which is ratio+1. 2856 2857 ! do iw=1,nw 2858 ! write(6,*) " freq, proc" 2859 ! write(6,*) iw, procb(iw,1) 2860 ! enddo 2861 endif 2862 2863 ! do iw=1,nw 2864 ! do iproc=0,nproc-1 2865 ! write(6,*) " freq, procs,?" 2866 ! write(6,*) iw,iproc,proct(iw,iproc) 2867 ! enddo 2868 ! enddo 2869 2870 2871 2872 end subroutine distrib_paral
m_green/fermi_green [ Functions ]
[ Top ] [ m_green ] [ Functions ]
NAME
fermi_green
FUNCTION
Compute Fermi level for DMFT or DFT.
COPYRIGHT
Copyright (C) 2006-2024 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 DFT fermi level max_iter : max number of iterations.
OUTPUT
fermie: output value
SOURCE
2985 subroutine fermi_green(cryst_struc,green,paw_dmft,pawang,self) 2986 2987 !Arguments ------------------------------------ 2988 !scalars 2989 type(crystal_t),intent(in) :: cryst_struc 2990 type(green_type),intent(inout) :: green 2991 type(paw_dmft_type), intent(inout) :: paw_dmft 2992 !type(MPI_type), intent(in) :: mpi_enreg 2993 type(pawang_type),intent(in) :: pawang 2994 type(self_type), intent(inout) :: self 2995 2996 !Local variables------------------------------- 2997 integer :: ierr_hh,opt_noninter,max_iter 2998 real(dp):: x_precision,f_precision,fermi_old 2999 ! real(dp) :: hx 3000 character(len=500) :: message 3001 !************************************************************************ 3002 ! 3003 write(message,'(a,8x,a)') ch10," == Compute Fermi level" 3004 call wrtout(std_out,message,'COLL') 3005 3006 !============= 3007 !headers 3008 !============= 3009 write(message,'(2a)') ch10, " |---Newton method to search Fermi level ------------|" 3010 call wrtout(std_out,message,'COLL') 3011 write(message,'(2a,f13.6)') ch10, " |--- Initial value for Fermi level",paw_dmft%fermie 3012 call wrtout(std_out,message,'COLL') 3013 3014 !======================================== 3015 !Define precision and nb of iterations 3016 !========================================= 3017 fermi_old=paw_dmft%fermie 3018 ierr_hh=0 3019 f_precision=paw_dmft%dmft_charge_prec 3020 !f_precision=0.01 3021 x_precision=tol5 3022 !if(option==1) then 3023 !f_precision=(erreursurlacharge)/100d0 3024 !else 3025 !f_precision=tol11 3026 !endif 3027 max_iter=1 ! for tests only 3028 !write(6,*) "for tests max_iter=1" 3029 max_iter=50 3030 opt_noninter=4 3031 3032 !===================== 3033 !Call newton method 3034 !===================== 3035 write(message,'(a,4x,a,e13.6)') ch10," Precision required :",f_precision 3036 call wrtout(std_out,message,'COLL') 3037 if (f_precision<10_dp) then 3038 call newton(cryst_struc,green,paw_dmft,pawang,self,& 3039 & paw_dmft%fermie,x_precision,max_iter,f_precision,ierr_hh,opt_noninter) 3040 end if 3041 3042 !=========================== 3043 !Deals with errors signals 3044 !=========================== 3045 if(ierr_hh==-314) then 3046 write(message,'(a)') "Warning, check Fermi level" 3047 call wrtout(std_out,message,'COLL') 3048 ! call abi_abort('COLL') 3049 write(message,'(2a,f13.6)') ch10, " |--- Final value for Fermi level (check)",paw_dmft%fermie 3050 call wrtout(std_out,message,'COLL') 3051 else if (ierr_hh==-123) then 3052 write(message,'(a,f13.6)') " Fermi level is put to",fermi_old 3053 paw_dmft%fermie=fermi_old 3054 call wrtout(std_out,message,'COLL') 3055 3056 ! ===================================== 3057 ! If fermi level search was successful 3058 ! ===================================== 3059 else 3060 write(message,'(a,4x,a,e13.6)') ch10," Precision achieved on Fermi Level :",x_precision 3061 call wrtout(std_out,message,'COLL') 3062 write(message,'(4x,a,e13.6)') " Precision achieved on number of electrons :",f_precision 3063 call wrtout(std_out,message,'COLL') 3064 write(message,'(2a,f13.6)') ch10, " |--- Final value for Fermi level",paw_dmft%fermie 3065 call wrtout(std_out,message,'COLL') 3066 end if 3067 3068 !======================================================== 3069 !Check convergence of fermi level during DMFT iterations 3070 !======================================================== 3071 if(paw_dmft%idmftloop>=2) then 3072 if(abs(paw_dmft%fermie-fermi_old).le.paw_dmft%dmft_fermi_prec) then 3073 ! write(message,'(a,8x,a,e9.2,a,8x,a,e12.5)') ch10,"|fermie(n)-fermie(n-1)|=<",paw_dmft%dmft_fermi_prec,ch10,& 3074 write(message,'(a,8x,a,e9.2,a,e9.2,a,8x,a,e12.5)') ch10,"|fermie(n)-fermie(n-1)|=",& 3075 & abs(paw_dmft%fermie-fermi_old),"<",paw_dmft%dmft_fermi_prec,ch10,& 3076 & "=> DMFT Loop: Fermi level is converged to:",paw_dmft%fermie 3077 call wrtout(std_out,message,'COLL') 3078 green%ifermie_cv=1 3079 else 3080 write(message,'(a,8x,a,2f12.5)') ch10,"DMFT Loop: Fermi level is not converged:",& 3081 & paw_dmft%fermie 3082 call wrtout(std_out,message,'COLL') 3083 green%ifermie_cv=0 3084 end if 3085 end if 3086 write(message,'(2a)') ch10, " |---------------------------------------------------|" 3087 call wrtout(std_out,message,'COLL') 3088 ! 3089 3090 !========================================================== 3091 !Recompute full green function including non diag elements 3092 !========================================================== 3093 call compute_green(cryst_struc,green,paw_dmft,pawang,0,self,opt_self=1,opt_nonxsum=1) 3094 call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=0,opt_ksloc=3) !,opt_nonxsum=1) 3095 3096 return 3097 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-2024 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
SOURCE
2499 subroutine fourier_fct(fw,ft,ldiag,ltau,opt_four,paw_dmft) 2500 2501 !Arguments ------------------------------------ 2502 !type 2503 logical,intent(in) :: ldiag 2504 integer,intent(in) :: ltau,opt_four 2505 type(paw_dmft_type), intent(in) :: paw_dmft 2506 complex(dpc), intent(inout) :: fw(paw_dmft%dmft_nwlo) 2507 complex(dpc), intent(inout) :: ft(ltau) 2508 2509 !local variables------------------------------- 2510 complex(dpc), allocatable :: splined_li(:) 2511 complex(dpc), allocatable :: tospline_li(:) 2512 ! complex(dpc), allocatable :: fw1(:) 2513 real(dp), allocatable :: ftr(:) 2514 real(dp) :: beta 2515 complex(dpc) :: xsto 2516 integer :: iflag,ifreq,itau,iwarn,log_direct 2517 character(len=500) :: message 2518 real(dp), allocatable :: omega_li(:) 2519 ! ********************************************************************* 2520 beta=one/paw_dmft%temp 2521 iflag=0 2522 log_direct=1 2523 2524 if(ldiag) iflag=1 2525 2526 ! == inverse fourier transform 2527 if(opt_four==-1) then 2528 2529 ABI_MALLOC(splined_li,(paw_dmft%dmft_nwli)) 2530 ! allocate(fw1(0:paw_dmft%dmft_nwlo-1)) 2531 if(paw_dmft%dmft_log_freq==1) then 2532 ABI_MALLOC(omega_li,(1:paw_dmft%dmft_nwli)) 2533 call construct_nwli_dmft(paw_dmft,paw_dmft%dmft_nwli,omega_li) 2534 call spline_c(paw_dmft%dmft_nwlo,paw_dmft%dmft_nwli,paw_dmft%omega_lo,& 2535 & omega_li,splined_li,fw) 2536 ABI_FREE(omega_li) 2537 else 2538 splined_li=fw 2539 endif 2540 call invfourier(splined_li,ft,paw_dmft%dmft_nwli,ltau,iflag,beta) 2541 ! deallocate(fw1) 2542 ABI_FREE(splined_li) 2543 2544 ! == direct fourier transform 2545 else if(opt_four==1) then 2546 2547 ABI_MALLOC(ftr,(ltau)) 2548 2549 iwarn=0 2550 do itau=1,ltau 2551 if(abs(aimag(ft(itau)))>tol12) then 2552 if(ldiag) then 2553 write(message,'(a,a,2(e15.4))') ch10,& 2554 & "green function is not real in imaginary time space",ft(itau) 2555 ABI_ERROR(message) 2556 else 2557 iwarn=iwarn+1 2558 ftr(itau)=real(ft(itau)) 2559 xsto=ft(itau) 2560 endif 2561 else 2562 ftr(itau)=real(ft(itau)) 2563 endif 2564 ! write(std_out,*) itau,ftr(itau) 2565 enddo 2566 2567 ABI_MALLOC(tospline_li,(paw_dmft%dmft_nwli)) 2568 if (log_direct==1) then 2569 ! do not have a physical meaning..because the log frequency is not 2570 ! one of the linear frequency. 2571 ! if log frequencies are also one of linear frequency, it should be good 2572 do ifreq=1, paw_dmft%dmft_nwlo 2573 call nfourier2(ftr,fw(ifreq),iflag,paw_dmft%omega_lo(ifreq),ltau,beta) 2574 enddo 2575 else 2576 call nfourier(ftr,tospline_li,iflag,paw_dmft%dmft_nwli-1,ltau,beta) 2577 do ifreq=1,paw_dmft%dmft_nwli 2578 write(1112,*) paw_dmft%temp*pi*real(2*ifreq-1,kind=dp),real(tospline_li(ifreq),kind=dp),aimag(tospline_li(ifreq)) 2579 !write(1112,*) paw_dmft%omega_li(ifreq),real(tospline_li(ifreq)),aimag(tospline_li(ifreq)) 2580 enddo 2581 if(paw_dmft%dmft_log_freq==1) then 2582 ABI_MALLOC(omega_li,(1:paw_dmft%dmft_nwli)) 2583 call construct_nwli_dmft(paw_dmft,paw_dmft%dmft_nwli,omega_li) 2584 call spline_c(paw_dmft%dmft_nwli,paw_dmft%dmft_nwlo,omega_li,& 2585 & paw_dmft%omega_lo,fw,tospline_li) 2586 ABI_FREE(omega_li) 2587 else 2588 fw=tospline_li 2589 endif 2590 endif 2591 2592 ABI_FREE(tospline_li) 2593 2594 ABI_FREE(ftr) 2595 if(iwarn>0) then 2596 write(message,'(a,a,2(e15.4))') ch10,& 2597 & "WARNING: off-diag green function is not real in imaginary time space",xsto 2598 call wrtout(std_out,message,'COLL') 2599 endif 2600 2601 endif 2602 2603 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
SOURCE
1869 subroutine fourier_green(cryst_struc,green,paw_dmft,pawang,opt_ksloc,opt_tw) 1870 1871 !Arguments ------------------------------------ 1872 !type 1873 type(crystal_t),intent(in) :: cryst_struc 1874 type(green_type),intent(inout) :: green 1875 !type(MPI_type), intent(in) :: mpi_enreg 1876 type(pawang_type),intent(in) :: pawang 1877 type(paw_dmft_type), intent(in) :: paw_dmft 1878 integer,intent(in) :: opt_ksloc,opt_tw ! fourier on ks or local 1879 1880 !local variables------------------------------- 1881 integer :: iatom,ib,ib1,ierr,ifreq,ikpt,im,im1,iparal,is,ispinor,ispinor1,itau 1882 integer :: mband,mbandc,myproc,natom,ndim,nkpt,nproc,nspinor,nsppol,spacecomm!,opt_four 1883 character(len=500) :: message 1884 ! complex(dpc):: ybcbeg,ybcend 1885 ! arrays 1886 complex(dpc), allocatable :: fw(:) 1887 complex(dpc), allocatable :: ft(:) 1888 type(green_type) :: green_temp 1889 ! ********************************************************************* 1890 ! ybcbeg=czero 1891 ! ybcend=czero 1892 1893 ! define spaceComm, myproc, and nproc from world communicator 1894 ! and mpi_enreg 1895 spacecomm=paw_dmft%spacecomm 1896 myproc=paw_dmft%myproc 1897 nproc=paw_dmft%nproc 1898 1899 ! Initialise integers 1900 mband = paw_dmft%mband 1901 mbandc = paw_dmft%mbandc 1902 nkpt = paw_dmft%nkpt 1903 nspinor = paw_dmft%nspinor 1904 nsppol = paw_dmft%nsppol 1905 natom = cryst_struc%natom 1906 1907 ! Only imaginary frequencies here 1908 if(green%w_type=="real") then 1909 message = 'fourier_green not implemented for real frequency' 1910 ABI_BUG(message) 1911 endif 1912 1913 ! Initialise temporary green function 1914 call init_green(green_temp,paw_dmft) 1915 call init_green_tau(green_temp,paw_dmft) 1916 1917 !green%oper(:)%matlu(1)%mat(1,1,1,1,1) 1918 ABI_MALLOC(fw,(green%nw)) 1919 ABI_MALLOC(ft,(green%dmftqmc_l)) 1920 1921 !============================================== 1922 ! == Inverse fourier transformation =========== 1923 !============================================== 1924 1925 if(opt_tw==-1) then 1926 1927 ! = For Kohn Sham green function 1928 !============================================== 1929 if(opt_ksloc==1.and.green%use_oper_tau_ks==1) then 1930 1931 do is = 1 , nsppol 1932 do ikpt = 1, nkpt 1933 do ib = 1, mbandc 1934 do ib1 = 1, mbandc 1935 do ifreq=1,green%nw 1936 fw(ifreq)=green%oper(ifreq)%ks(is,ikpt,ib,ib1) 1937 enddo 1938 call fourier_fct(fw,ft,ib==ib1,green%dmftqmc_l,-1,paw_dmft) ! inverse fourier 1939 do itau=1,green%dmftqmc_l 1940 green_temp%oper_tau(itau)%ks(is,ikpt,ib,ib1)=ft(itau) 1941 enddo 1942 if(ib==ib1) then 1943 green%occup_tau%ks(is,ikpt,ib,ib1)=ft(1)+one 1944 else 1945 green%occup_tau%ks(is,ikpt,ib,ib1)=ft(1) 1946 endif 1947 enddo ! ib1 1948 enddo ! ib 1949 enddo ! ikpt 1950 enddo ! isppol 1951 1952 ! = Post-treatment: necessary in the case of nspinor==2, but valid anywhere 1953 ! because G(tau)=G_{nu,nu'}(tau)+[G_{nu,nu'}(tau)]* 1954 1955 do is = 1 , nsppol 1956 do ikpt = 1, nkpt 1957 do ib = 1, mbandc 1958 do ib1 = 1, mbandc 1959 do itau=1,green%dmftqmc_l 1960 green%oper_tau(itau)%ks(is,ikpt,ib,ib1)=& 1961 & (green_temp%oper_tau(itau)%ks(is,ikpt,ib,ib1)+ & 1962 & conjg(green_temp%oper_tau(itau)%ks(is,ikpt,ib1,ib)))/two 1963 if(ib==ib1) then 1964 green%occup_tau%ks(is,ikpt,ib,ib1)=green%oper_tau(1)%ks(is,ikpt,ib,ib1)+one 1965 else 1966 green%occup_tau%ks(is,ikpt,ib,ib1)=green%oper_tau(1)%ks(is,ikpt,ib,ib1) 1967 endif 1968 enddo 1969 enddo ! ib1 1970 enddo ! ib 1971 enddo ! ikpt 1972 enddo ! isppol 1973 call loc_oper(green%occup_tau,paw_dmft,1) 1974 call sym_matlu(cryst_struc,green%occup_tau%matlu,pawang,paw_dmft) 1975 write(message,'(a,a,i10,a)') ch10," green%occup_tau%matlu from green_occup_tau%ks" 1976 call wrtout(std_out,message,'COLL') 1977 call print_matlu(green%occup_tau%matlu,natom,prtopt=3) 1978 endif 1979 1980 ! = For local green function 1981 !============================================== 1982 if(opt_ksloc ==2) then 1983 1984 iparal=0 1985 do iatom=1, natom 1986 if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then 1987 ndim=2*green%oper(1)%matlu(iatom)%lpawu+1 1988 do is = 1 , nsppol 1989 do ispinor = 1, nspinor 1990 do ispinor1 = 1, nspinor 1991 do im=1,ndim 1992 do im1=1,ndim 1993 iparal=iparal+1 1994 if(mod(iparal-1,nproc)==myproc) then 1995 do ifreq=1,green%nw 1996 fw(ifreq)=green%oper(ifreq)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1) 1997 enddo 1998 ! inverse fourier 1999 ! write(std_out,'(a)') "fourierbeforeposttreatement,ispinor,ispinor1,is,im,im1" 2000 ! write(std_out,'(a,5i4,f12.5,f12.5)') "fourier",ispinor,ispinor1,is,im,im1 2001 ! write(std_out,'(a,e12.5,e12.5)')& 2002 ! &"green%oper(4)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)"& 2003 !& ,green%oper(4)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1) 2004 call fourier_fct(fw,ft,(im==im1).and.(ispinor==ispinor1),green%dmftqmc_l,-1,paw_dmft) 2005 do itau=1,green%dmftqmc_l 2006 green_temp%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=ft(itau) 2007 ! write(std_out,*) itau,green_temp%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1) 2008 enddo 2009 ! if((im==im1).and.(ispinor==ispinor1)) then 2010 ! green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=ft(1)+one 2011 ! else 2012 ! green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=ft(1) 2013 ! endif 2014 endif ! iparal 2015 enddo 2016 enddo 2017 enddo ! ispinor1 2018 enddo ! ispinor 2019 enddo ! isppol 2020 endif ! lpawu.ne.-1 2021 enddo ! iatom 2022 2023 ! Parallelisation must be finished here, because in the post 2024 ! treatment, value from different proc will be mixed. 2025 call xmpi_barrier(spacecomm) 2026 do iatom=1,natom 2027 do itau=1,green%dmftqmc_l 2028 call xmpi_sum(green_temp%oper_tau(itau)%matlu(iatom)%mat,spacecomm,ierr) 2029 enddo 2030 enddo 2031 call xmpi_barrier(spacecomm) 2032 2033 ! = Post-treatment: necessary in the case of nspinor==2, but valid anywhere 2034 ! because G(tau)=G_{LL'}^{sigma,sigma'}(tau)+[G_{L'L}^{sigma',sigma}(tau)]* 2035 2036 if(nspinor>0) Then 2037 do iatom=1, natom 2038 if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then 2039 ndim=2*green%oper(1)%matlu(iatom)%lpawu+1 2040 do is = 1 , nsppol 2041 do ispinor = 1, nspinor 2042 do ispinor1 = 1, nspinor 2043 do im=1,ndim 2044 do im1=1,ndim 2045 ! write(std_out,'(a,5i4,f12.5,f12.5)') "fourier -1",ispinor,ispinor1,is,im,im1 2046 do itau=1,green%dmftqmc_l 2047 green%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=& 2048 & ((green_temp%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)+ & 2049 & conjg(green_temp%oper_tau(itau)%matlu(iatom)%mat(im1,im,is,ispinor1,ispinor))))/two 2050 ! write(std_out,*) itau,green%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1) 2051 if((im==im1).and.(ispinor==ispinor1)) then 2052 green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=& 2053 & green%oper_tau(1)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)+one 2054 else 2055 green%occup_tau%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=& 2056 & green%oper_tau(1)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1) 2057 endif ! test diag 2058 enddo ! itau 2059 enddo ! im1 2060 enddo ! im 2061 enddo ! ispinor1 2062 enddo ! ispinor 2063 enddo ! isppol 2064 endif 2065 enddo ! iatom 2066 endif ! nspinor>0 2067 endif ! opt_ksloc=2 2068 endif ! opt_tw==-1 2069 2070 2071 !============================================== 2072 ! == Direct fourier transformation 2073 !============================================== 2074 ! todo_ba ft useful only for diagonal elements ... 2075 2076 if(opt_tw==1) then 2077 2078 ! = For local green function 2079 !============================================== 2080 if(opt_ksloc ==2) then 2081 2082 iparal=0 2083 do iatom=1, natom 2084 ! put to zero (usefull at least for parallelism) 2085 do ifreq=1,green%nw 2086 green%oper(ifreq)%matlu(iatom)%mat=czero 2087 enddo 2088 if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then 2089 ndim=2*green%oper(1)%matlu(iatom)%lpawu+1 2090 do is = 1 , nsppol 2091 do ispinor = 1, nspinor 2092 do ispinor1 = 1, nspinor 2093 do im=1,ndim 2094 do im1=1,ndim 2095 iparal=iparal+1 2096 if(mod(iparal-1,nproc)==myproc) then 2097 do itau=1,green%dmftqmc_l 2098 ft(itau)=green%oper_tau(itau)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1) 2099 enddo 2100 ! if(im==1.and.im1==1) write(std_out,*) "ft(itau=1)",is,ft(1) ! debug 2101 call fourier_fct(fw,ft,(im==im1).and.(ispinor==ispinor1),green%dmftqmc_l,1,paw_dmft) 2102 do ifreq=1,green%nw 2103 green%oper(ifreq)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=fw(ifreq) 2104 enddo 2105 ! if(im==1.and.im1==1) write(std_out,*) "fw(ifreq=1)",is,fw(1) ! debug 2106 endif 2107 enddo 2108 enddo 2109 enddo ! ispinor1 2110 enddo ! ispinor 2111 enddo ! isppol 2112 endif ! lpawu=/-1 2113 enddo ! iatom 2114 call xmpi_barrier(spacecomm) 2115 do iatom=1,natom 2116 do ifreq=1,green%nw 2117 call xmpi_sum(green%oper(ifreq)%matlu(iatom)%mat,spacecomm,ierr) 2118 enddo 2119 enddo 2120 endif ! opt_ksloc=2 2121 endif ! opt_tw==-1 2122 ABI_FREE(fw) 2123 ABI_FREE(ft) 2124 call destroy_green_tau(green_temp) 2125 call destroy_green(green_temp) 2126 2127 end subroutine fourier_green
m_green/green_type [ Types ]
NAME
green_type
FUNCTION
This structured datatype contains the necessary data
SOURCE
83 type, public :: green_type ! for each atom 84 85 integer :: dmft_nwlo 86 ! dmft frequencies 87 88 character(len=4) :: w_type 89 ! type of frequencies used 90 91 character(len=12) :: whichgreen 92 ! describe the type of green function computed (DFT, DMFT, ..) 93 94 integer :: nw 95 ! dmft frequencies 96 97 integer :: fileprt_w 98 ! 1 if file created 99 100 integer :: fileprt_tau 101 ! 1 if file created 102 103 integer :: dmft_nwli 104 ! dmft frequencies 105 106 integer :: dmftqmc_l 107 ! number of time slices for QMC 108 109 integer :: ifermie_cv 110 111 integer :: ichargeloc_cv 112 113 integer :: use_oper_tau_ks 114 ! 0 do not use oper_tau_ks 115 ! 1 use oper_tau_ks 116 117 real(dp) :: charge_ks 118 ! Total charge computed from ks orbitals 119 120 integer :: has_charge_matlu_solver 121 ! =0 charge_matlu_solver not allocated 122 ! =1 charge_matlu_solver is allocated 123 ! =2 charge_matlu_solver is calculated (ie calculation of LOCAL CORRELATED occupations is done from 124 ! solver green function) 125 126 integer :: has_greenmatlu_xsum 127 ! =1 green%oper%matlu xsum in compute_green 128 ! =0 green%oper%matlu non xsumed in compute_green 129 ! used in integrate_green to checked that green function was computed in 130 ! integrate_green. 131 132 integer :: has_charge_matlu 133 ! =2 if calculation of LOCAL CORRELATED occupations is done from 134 135 integer :: has_charge_matlu_prev 136 ! =0 charge_matlu_prev not allocated 137 ! =1 charge_matlu_prev is allocated 138 ! =2 charge_matlu_prev is calculated (ie calculation of LOCAL CORRELATED occupations is done from 139 ! solver green function) 140 141 integer, allocatable :: procb(:,:) 142 143 integer, allocatable :: proct(:,:) 144 145 real(dp), allocatable :: charge_matlu_prev(:,:) 146 ! Total charge on correlated orbitals from previous iteration 147 148 real(dp), allocatable :: charge_matlu(:,:) 149 ! Total charge on correlated orbitals 150 ! todo_ba name of charge_matlu is misleading: should be changed 151 152 real(dp), allocatable :: charge_matlu_solver(:,:) 153 ! Total charge on correlated orbitals obtained from solver by 154 ! integration over frequencies. 155 156 real(dp), allocatable :: tau(:) 157 ! value of time in imaginary space 158 159 real(dp), pointer :: omega(:) => null() 160 ! value of frequencies 161 162 real(dp), allocatable :: ecorr_qmc(:) 163 ! Correlation energy for a given atom in qmc 164 165 type(oper_type), allocatable :: oper(:) 166 ! green function in different basis 167 168 type(oper_type), allocatable :: oper_tau(:) 169 ! green function in different basis 170 171 type(oper_type) :: occup 172 ! occupation in different basis 173 174 type(oper_type) :: occup_tau 175 ! occupation in different basis 176 177 end type green_type 178 179 !---------------------------------------------------------------------- 180 181 182 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=information 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
SOURCE
1778 subroutine icip_green(char1,cryst_struc,green,paw_dmft,pawang,pawprtvol,self,opt_self) 1779 1780 !Arguments ------------------------------------ 1781 !type 1782 type(crystal_t),intent(in) :: cryst_struc 1783 !type(MPI_type), intent(in) :: mpi_enreg 1784 type(green_type),intent(inout) :: green 1785 type(pawang_type),intent(in) :: pawang 1786 type(paw_dmft_type), intent(inout) :: paw_dmft 1787 type(self_type), intent(inout) :: self 1788 integer, intent(in) :: pawprtvol 1789 character(len=*), intent(in) :: char1 1790 integer, optional, intent(in) :: opt_self 1791 !Local variables ------------------------------------ 1792 integer :: option,optself,optlocks,prtopt_for_integrate_green,opt_nonxsum_icip 1793 character(len=500) :: message 1794 ! ********************************************************************* 1795 green%whichgreen="DFT" 1796 prtopt_for_integrate_green=2 1797 1798 if(present(opt_self)) then 1799 optself=opt_self 1800 else 1801 optself=0 1802 endif 1803 opt_nonxsum_icip=1 1804 if(paw_dmft%dmftcheck==1) then ! for fourier_green 1805 ! opt_nonxsum_icip=0 1806 endif 1807 1808 call init_green(green,paw_dmft) 1809 1810 ! useless test ok 1811 ! call printocc_green(green,1,paw_dmft,2) 1812 ! write(std_out,*)" printocc_green zero finished " 1813 1814 !== Compute green%oper(:)%ks 1815 !== Deduce green%oper(:)%matlu(:)%mat 1816 call compute_green(cryst_struc,green,paw_dmft,& 1817 & pawang,pawprtvol,self,optself,opt_nonxsum=opt_nonxsum_icip) 1818 if(paw_dmft%dmft_prgn>=1.and.paw_dmft%dmft_prgn<=2) then 1819 optlocks=paw_dmft%dmft_prgn*2+1 ! if dmft_prgn==2 => do not print 1820 if(paw_dmft%lpsichiortho==1.and.pawprtvol>-100) then 1821 call print_green(char1,green,optlocks,paw_dmft,pawprtvol) 1822 ! call print_green('inicip',green,1,paw_dmft,pawprtvol=1,opt_wt=1) 1823 endif 1824 endif 1825 1826 !== Integrate green%oper(:)%ks 1827 !== Integrate green%oper(:)%matlu(:)%mat 1828 call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt_for_integrate_green,opt_ksloc=3)!,opt_nonxsum=opt_nonxsum_icip) 1829 !== Print green%oper(:)%ks 1830 !== Print green%oper(:)%matlu(:)%mat 1831 if(char1=="DFT") then 1832 option=1 1833 if(self%oper(1)%matlu(1)%lpawu/=-1) then 1834 if(abs(real(self%oper(1)%matlu(1)%mat(1,1,1,1,1)))>tol7) then 1835 ! todo_ab: generalise this 1836 write(message,'(a,a,2(e15.4))') ch10,& 1837 & "Warning: a DFT calculation is carried out and self is not zero" 1838 call wrtout(std_out,message,'COLL') 1839 ! call abi_abort('COLL') 1840 endif 1841 endif 1842 else 1843 option=5 1844 endif 1845 if(pawprtvol>-100) then 1846 call printocc_green(green,option,paw_dmft,3,chtype=char1) 1847 end if 1848 1849 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
SOURCE
204 subroutine init_green(green,paw_dmft,opt_oper_ksloc,wtype) 205 206 !Arguments ------------------------------------ 207 !scalars 208 !type 209 type(green_type), intent(out) :: green 210 type(paw_dmft_type), intent(in) :: paw_dmft 211 integer, optional, intent(in) :: opt_oper_ksloc 212 character(len=4), optional :: wtype 213 !local variables ------------------------------------ 214 integer :: ifreq,nw,optoper_ksloc,nw_perproc 215 !************************************************************************ 216 if(present(opt_oper_ksloc)) then 217 optoper_ksloc=opt_oper_ksloc 218 else 219 optoper_ksloc=2 220 endif 221 if(present(wtype)) then 222 green%w_type=wtype 223 else 224 green%w_type="imag" 225 endif 226 227 228 if(green%w_type=="imag") then 229 nw=paw_dmft%dmft_nwlo 230 green%omega=>paw_dmft%omega_lo 231 else if(green%w_type=="real") then 232 nw=size(paw_dmft%omega_r) 233 green%omega=>paw_dmft%omega_r 234 endif 235 236 green%dmft_nwlo=paw_dmft%dmft_nwlo 237 green%dmft_nwli=paw_dmft%dmft_nwli 238 green%charge_ks=zero 239 green%has_charge_matlu=0 240 ABI_MALLOC(green%charge_matlu,(paw_dmft%natom,paw_dmft%nsppol+1)) 241 green%charge_matlu=zero 242 green%has_charge_matlu=1 243 green%has_greenmatlu_xsum=0 244 245 green%has_charge_matlu_solver=0 246 ABI_MALLOC(green%charge_matlu_solver,(paw_dmft%natom,paw_dmft%nsppol+1)) 247 green%charge_matlu_solver=zero 248 green%has_charge_matlu_solver=1 249 250 green%has_charge_matlu_prev=0 251 ABI_MALLOC(green%charge_matlu_prev,(paw_dmft%natom,paw_dmft%nsppol+1)) 252 green%charge_matlu_prev=zero 253 green%has_charge_matlu_prev=1 254 255 call init_oper(paw_dmft,green%occup,opt_ksloc=3) 256 257 ! built simple arrays to distribute the tasks in compute_green. 258 ABI_MALLOC(green%procb,(nw,paw_dmft%nkpt)) 259 ABI_MALLOC(green%proct,(nw,0:paw_dmft%nproc-1)) 260 261 call distrib_paral(paw_dmft%nkpt,paw_dmft%nproc,nw,nw_perproc,green%procb,green%proct) 262 green%nw=nw 263 264 ! need to distribute memory over frequencies 265 266 !! begin of temporary modificatios 267 ! ABI_MALLOC(green%oper,(green%nw_perproc)) 268 ! do ifreq=1,green%nw_perproc 269 ! call init_oper(paw_dmft,green%oper(ifreq),opt_ksloc=optoper_ksloc) 270 ! enddo 271 ! 272 ! do ifreq=1,green%nw 273 ! if(green%proct(ifreq,myproc)==1) then 274 ! do ikpt = 1 , paw_dmft%nkpt 275 ! if (green%procb(ifreq,ikpt)==myproc) then 276 ! endif 277 ! enddo ! ikpt 278 ! endif ! parallelisation 279 ! enddo ! ifreq 280 ! 281 ! 282 !! end of temporary modificatios 283 ABI_MALLOC(green%oper,(green%nw)) 284 do ifreq=1,green%nw 285 call init_oper(paw_dmft,green%oper(ifreq),opt_ksloc=optoper_ksloc) 286 enddo 287 green%ifermie_cv=0 288 green%ichargeloc_cv=0 289 290 if(paw_dmft%dmft_solv>=4) then 291 ABI_MALLOC(green%ecorr_qmc,(paw_dmft%natom)) 292 end if 293 if(paw_dmft%dmft_solv>=4) green%ecorr_qmc(:)=zero 294 295 296 green%fileprt_w=0 297 green%fileprt_tau=0 298 299 300 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
SOURCE
319 subroutine init_green_tau(green,paw_dmft,opt_ksloc) 320 321 !Arguments ------------------------------------ 322 !scalars 323 !type 324 type(green_type), intent(inout) :: green 325 type(paw_dmft_type), intent(in) :: paw_dmft 326 integer, optional, intent(in) :: opt_ksloc 327 !local variables ------------------------------------ 328 integer :: optksloc 329 integer :: itau 330 !************************************************************************ 331 if(present(opt_ksloc)) then 332 optksloc=opt_ksloc 333 else 334 optksloc=3 335 endif 336 green%use_oper_tau_ks=0 337 if(green%use_oper_tau_ks==0) then 338 optksloc=2 339 endif 340 341 green%dmftqmc_l=paw_dmft%dmftqmc_l 342 ABI_MALLOC(green%tau,(green%dmftqmc_l)) 343 do itau=1,green%dmftqmc_l 344 green%tau(itau)=float(itau-1)/float(green%dmftqmc_l)/paw_dmft%temp 345 enddo 346 347 call init_oper(paw_dmft,green%occup_tau,optksloc) 348 349 ABI_MALLOC(green%oper_tau,(paw_dmft%dmftqmc_l)) 350 do itau=1,green%dmftqmc_l 351 call init_oper(paw_dmft,green%oper_tau(itau),optksloc) 352 enddo 353 354 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-2024 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
SOURCE
2379 subroutine int_fct(ff,ldiag,option,paw_dmft,integral,procb,myproc) 2380 2381 !Arguments ------------------------------------ 2382 !type 2383 logical,intent(in) :: ldiag 2384 integer,intent(in) :: option 2385 complex(dpc),intent(out) :: integral 2386 type(paw_dmft_type), intent(in) :: paw_dmft 2387 complex(dpc), intent(in) :: ff(paw_dmft%dmft_nwlo) 2388 integer, optional, intent(in) :: procb(paw_dmft%dmft_nwlo) 2389 integer, optional, intent(in) :: myproc 2390 2391 !local variables------------------------------- 2392 logical, allocatable :: procb2(:) 2393 character(len=500) :: message 2394 integer :: ifreq 2395 ! ********************************************************************* 2396 ABI_MALLOC(procb2,(paw_dmft%dmft_nwlo)) 2397 if(present(procb).and.present(myproc)) then 2398 do ifreq=1,paw_dmft%dmft_nwlo 2399 procb2(ifreq)=(procb(ifreq)==myproc) 2400 enddo 2401 else if(present(procb).and..not.present(myproc)) then 2402 write(message,'(a,a,2(e15.4))') ch10,& 2403 & "BUG: procb is present and not myproc in int_fct" 2404 ABI_BUG(message) 2405 else if(.not.present(procb).and.present(myproc)) then 2406 write(message,'(a,a,2(e15.4))') ch10,& 2407 & "BUG: procb is not present and myproc is in int_fct" 2408 ABI_BUG(message) 2409 else 2410 do ifreq=1,paw_dmft%dmft_nwlo 2411 procb2(ifreq)=(1==1) 2412 enddo 2413 endif 2414 2415 integral=czero 2416 if(ldiag) then 2417 2418 if(option==1) then ! nspinor==1 2419 do ifreq=1,paw_dmft%dmft_nwlo 2420 if(procb2(ifreq)) then 2421 ! write(500,*) paw_dmft%omega_lo(ifreq),real(ff(ifreq)),imag(ff(ifreq)) 2422 integral=integral+2.d0*paw_dmft%temp * & 2423 & real( ff(ifreq)-one / ( j_dpc*paw_dmft%omega_lo(ifreq) ) ) * & 2424 & paw_dmft%wgt_wlo(ifreq) 2425 endif 2426 enddo 2427 if(procb2(paw_dmft%dmft_nwlo)) integral=integral+half 2428 ! integral=integral+half 2429 ! the if is here, to count only one time this correction 2430 endif 2431 2432 if(option==2) then ! nspinor==2 2433 do ifreq=1,paw_dmft%dmft_nwlo 2434 if(procb2(ifreq)) then 2435 integral=integral+2.d0*paw_dmft%temp * & 2436 & ( ff(ifreq)-one / ( j_dpc*paw_dmft%omega_lo(ifreq) ) ) * & 2437 & paw_dmft%wgt_wlo(ifreq) 2438 endif 2439 enddo 2440 if(procb2(paw_dmft%dmft_nwlo)) integral=integral+half 2441 ! integral=integral+half 2442 ! the if is here, to count only one time this correction 2443 endif 2444 2445 2446 else ! ldiag 2447 2448 ! write(std_out,*) "nondiag" 2449 if(option==1) then 2450 do ifreq=1,paw_dmft%dmft_nwlo 2451 if(procb2(ifreq)) then 2452 integral=integral+2.d0*paw_dmft%temp * & 2453 & real( ff(ifreq) ) * & 2454 & paw_dmft%wgt_wlo(ifreq) 2455 endif 2456 enddo 2457 endif 2458 if(option==2) then 2459 do ifreq=1,paw_dmft%dmft_nwlo 2460 if(procb2(ifreq)) then 2461 integral=integral+2.d0*paw_dmft%temp * & 2462 & ff(ifreq) * & 2463 & paw_dmft%wgt_wlo(ifreq) 2464 endif 2465 enddo 2466 endif 2467 endif ! ldiag 2468 ABI_FREE(procb2) 2469 2470 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
SOURCE
1326 subroutine integrate_green(cryst_struc,green,paw_dmft& 1327 & ,pawang,prtopt,opt_ksloc,opt_after_solver,opt_diff,opt_nonxsum) 1328 1329 !Arguments ------------------------------------ 1330 !type 1331 type(crystal_t),intent(in) :: cryst_struc 1332 type(green_type),intent(inout) :: green 1333 !type(MPI_type), intent(in) :: mpi_enreg 1334 type(pawang_type),intent(in) :: pawang 1335 type(paw_dmft_type), intent(inout) :: paw_dmft 1336 integer, intent(in) :: prtopt 1337 integer, optional, intent(in) :: opt_ksloc 1338 integer, optional, intent(in) :: opt_after_solver 1339 integer, optional, intent(in) :: opt_diff 1340 integer, optional, intent(in) :: opt_nonxsum 1341 1342 !local variables------------------------------- 1343 real(dp) :: tsec(2) 1344 integer :: iatom,ib,ib1,icomp_chloc,ifreq,ikpt,im,im1,ispinor,ispinor1,is 1345 integer :: mband,mbandc,myproc,natom,ndim,nkpt,nproc,nspinor 1346 integer :: nsppol,option 1347 integer :: optksloc,spacecomm,optaftsolv,optnonxsum 1348 complex(dpc) :: integral 1349 character(len=500) :: message 1350 complex(dpc), allocatable :: ff(:) 1351 type(matlu_type), allocatable :: matlu_temp(:) 1352 type(oper_type) :: occup_temp 1353 real(dp) :: diff_chloc 1354 ! real(dp), allocatable :: charge_loc_old(:,:) 1355 ! type(oper_type) :: oper_c 1356 ! ********************************************************************* 1357 1358 DBG_ENTER("COLL") 1359 call timab(625,1,tsec) 1360 1361 if(prtopt>0) then 1362 write(message,'(2a,i3,13x,a)') ch10,' === Integrate green function' 1363 call wrtout(std_out,message,'COLL') 1364 endif 1365 if(green%w_type=="real") then 1366 message = 'integrate_green not implemented for real frequency' 1367 ABI_BUG(message) 1368 endif 1369 if(present(opt_nonxsum)) then 1370 optnonxsum=opt_nonxsum 1371 else 1372 optnonxsum=0 1373 endif 1374 1375 ! Initialise spaceComm, myproc, and master 1376 spacecomm=paw_dmft%spacecomm 1377 myproc=paw_dmft%myproc 1378 nproc=paw_dmft%nproc 1379 1380 1381 ! Initialise integers 1382 mband = paw_dmft%mband 1383 mbandc = paw_dmft%mbandc 1384 nkpt = paw_dmft%nkpt 1385 nspinor = paw_dmft%nspinor 1386 nsppol = paw_dmft%nsppol 1387 natom = paw_dmft%natom 1388 1389 ! Initialize green%oper before calculation (important for xmpi_sum) 1390 ! allocate(charge_loc_old(paw_dmft%natom,paw_dmft%nsppol+1)) 1391 ! if(.not.present(opt_diff)) then ! if integrate_green is called in m_dmft after calculation of self 1392 ! charge_loc_old=green%charge_matlu 1393 ! endif 1394 icomp_chloc=0 1395 1396 ! Choose what to compute 1397 if(present(opt_ksloc)) then 1398 optksloc=opt_ksloc 1399 else 1400 optksloc=3 1401 endif 1402 if(present(opt_after_solver)) then 1403 optaftsolv=opt_after_solver 1404 else 1405 optaftsolv=0 1406 endif 1407 if(optaftsolv==1.and.abs(optksloc)/=2) then 1408 message = "integration of ks green function should not be done after call to solver : it has not been computed" 1409 ABI_BUG(message) 1410 endif 1411 if(abs(optksloc)>=2.and.green%has_greenmatlu_xsum==0) then 1412 write(message,'(4a)') ch10,& 1413 & "BUG: integrate_green is asked to integrate local green function",ch10,& 1414 & " and local green function was non broadcasted in compute_green" 1415 ABI_BUG(message) 1416 endif 1417 1418 ! Allocations 1419 ABI_MALLOC(matlu_temp,(natom)) 1420 call init_matlu(natom,nspinor,nsppol,green%oper(1)%matlu(:)%lpawu,matlu_temp) 1421 1422 ABI_MALLOC(ff,(green%nw)) 1423 1424 ! ================================================= 1425 ! == Integrate Local Green function =============== 1426 if(abs(optksloc)/2==1) then ! optksloc=2 or 3 1427 ! ================================================= 1428 call zero_matlu(green%occup%matlu,green%occup%natom) 1429 1430 ! == Calculation of \int{G_{LL'}{\sigma\sigma',s}(R)(i\omega_n)} 1431 if(paw_dmft%lpsichiortho==1) then 1432 ! - Calculation of frequency sum over positive frequency 1433 if (nspinor==1) option=1 1434 if (nspinor==2) option=2 1435 do iatom=1, natom 1436 ndim=2*green%oper(1)%matlu(iatom)%lpawu+1 1437 if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then 1438 do ispinor1 = 1, nspinor 1439 do ispinor = 1, nspinor 1440 do is = 1 , nsppol 1441 do im=1,ndim 1442 do im1=1,ndim 1443 do ifreq=1,green%nw 1444 ff(ifreq)= & 1445 & green%oper(ifreq)%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1) 1446 ! write(std_out,*) green%omega(ifreq),ff(ifreq)," integrate green fw_lo" 1447 !if(present(iii).and.im==1.and.im1==1) write(std_out_default,*) ch10,ifreq,ff(ifreq),"#ff" 1448 enddo 1449 ! call int_fct(ff,(im==im1).and.(ispinor==ispinor1),& 1450 !& option,paw_dmft,integral) 1451 call int_fct(ff,(im==im1).and.(ispinor==ispinor1),& 1452 & 2,paw_dmft,integral) ! test_1 1453 green%occup%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)=integral 1454 !if(present(iii).and.im==1.and.im1==1) write(std_out_default,*) ch10,'integral',im,im1,ifreq,integral 1455 ! if(im==2.and.im1==5.and.is==1.and.iatom==1) then 1456 ! write(std_out,*) " occup %matlu(1)%mat(2,5,1,1,1)",green%occup%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1) 1457 ! endif 1458 enddo 1459 enddo 1460 enddo ! ispinor1 1461 enddo ! ispinor 1462 enddo ! is 1463 matlu_temp(iatom)%mat=green%occup%matlu(iatom)%mat 1464 endif ! lpawu=/-1 1465 enddo ! iatom 1466 1467 ! Print density matrix if prtopt high 1468 if(abs(prtopt)>2) then 1469 write(message,'(a,a,i10,a)') ch10," = green%occup%matlu from int(gloc(w))" 1470 call wrtout(std_out,message,'COLL') 1471 call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=-1) 1472 endif 1473 1474 ! - Symetrise: continue sum over k-point: Full BZ 1475 call sym_matlu(cryst_struc,green%occup%matlu,pawang,paw_dmft) 1476 if(abs(prtopt)>2) then 1477 write(message,'(a,a,i10,a)') ch10,& 1478 & " = green%occup%matlu from int(gloc(w)) with symetrization" 1479 call wrtout(std_out,message,'COLL') 1480 call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=-1) 1481 endif 1482 call sym_matlu(cryst_struc,matlu_temp,pawang,paw_dmft) 1483 1484 ! - Post-treatment for summation over negative and positive frequencies: 1485 ! necessary in the case of nspinor==2 AND nspinor==1, but valid anywhere 1486 ! N(ll'sigmasigma')= (N(ll'sigmasigma')+ N*(l'lsigma'sigma))/2 1487 ! because [G_{LL'}^{sigma,sigma'}(iomega_n)]*= G_{L'L}^{sigma',sigma}(-iomega_n) 1488 if(nspinor>=1) Then 1489 do iatom=1, natom 1490 ndim=2*green%oper(1)%matlu(iatom)%lpawu+1 1491 if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then 1492 do ispinor1 = 1, nspinor 1493 do ispinor = 1, nspinor 1494 do is = 1 , nsppol 1495 do im=1,ndim 1496 do im1=1,ndim 1497 green%occup%matlu(iatom)%mat(im,im1,is,ispinor,ispinor1)= & 1498 & (matlu_temp(iatom)%mat(im,im1,is,ispinor,ispinor1)+ & 1499 & conjg(matlu_temp(iatom)%mat(im1,im,is,ispinor1,ispinor)))/two 1500 enddo 1501 enddo 1502 enddo ! ispinor1 1503 enddo ! ispinor 1504 enddo ! isppol 1505 matlu_temp(iatom)%mat=green%occup%matlu(iatom)%mat 1506 endif 1507 enddo ! iatom 1508 endif 1509 if(abs(prtopt)>2) then 1510 write(message,'(a,a,i10,a)') ch10,& 1511 & " = green%occup%matlu from int(gloc(w)) symetrized with post-treatment" 1512 call wrtout(std_out,message,'COLL') 1513 call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=-1) 1514 endif 1515 if(optaftsolv==0) then 1516 call trace_oper(green%occup,green%charge_ks,green%charge_matlu,2) 1517 green%has_charge_matlu=2 1518 green%has_charge_matlu_solver=0 1519 icomp_chloc=1 1520 else if(optaftsolv==1.and.paw_dmft%dmft_solv/=4) then 1521 ! else if(optaftsolv==1) then 1522 if(paw_dmft%dmft_solv==4) then 1523 write(message,'(a,a,a)') ch10,& 1524 & " = WARNING: Double Counting will be computed with solver charge.." ,& 1525 & "might be problematic with Hirsch Fye QMC" 1526 call wrtout(std_out,message,'COLL') 1527 endif 1528 ! This is done only when called from impurity_solver with solver 1529 ! green function. (And if QMC is NOT used). 1530 if(paw_dmft%dmft_solv>=4.and.green%has_charge_matlu_solver/=2) then 1531 write(message,'(a,a,i3)') ch10,& 1532 & " = BUG : has_charge_matlu_solver should be 2 and is",& 1533 & green%has_charge_matlu_solver 1534 ABI_BUG(message) 1535 endif 1536 if(paw_dmft%dmft_solv<=4) then 1537 call trace_oper(green%occup,green%charge_ks,green%charge_matlu_solver,2) 1538 green%has_charge_matlu_solver=2 1539 endif 1540 endif 1541 else 1542 write(message,'(a,4x,a,a,a,4x,a)') ch10,& 1543 & " Local basis is not (yet) orthonormal:",& 1544 & " local green function is thus not integrated",ch10,& 1545 & " Local occupations are computed from KS occupations" 1546 call wrtout(std_out,message,'COLL') 1547 endif 1548 1549 endif ! optksloc 1550 ! ================================================= 1551 1552 1553 1554 ! ================================================= 1555 ! == Integrate Kohn Sham Green function =========== 1556 if(mod(abs(optksloc),2)==1) then ! optksloc=1 or 3 or -1 1557 ! green%occup%ks=czero ! important for xmpi_sum 1558 !! ================================================= 1559 !! == Calculation of \int{G_{\nu\nu'}{k,s}(i\omega_n)} 1560 call init_oper(paw_dmft,occup_temp) 1561 ! ff=czero 1562 ! do is = 1 , nsppol 1563 ! do ikpt = 1, nkpt 1564 !!020212 if(mod(ikpt-1,nproc)==myproc) then 1565 !!! print'(A,3I4)', "P ",myproc,is,ikpt 1566 ! do ib = 1, mbandc 1567 ! do ib1 = 1, mbandc 1568 ! if(optksloc==-1.and.(ib/=ib1)) cycle 1569 ! ff(:)=czero 1570 ! do ifreq=1,green%nw 1571 ! !! the following line is a quick and dirty tricks and should be removed when the data will 1572 ! !! be correctly distributed. 1573 !!! if((optnonxsum==1).and.(green%procb(ifreq,ikpt)==myproc)).or.(optnonxsum==0)) then 1574 ! if(green%procb(ifreq,ikpt)==myproc) then 1575 ! ff(ifreq)=green%oper(ifreq)%ks(is,ikpt,ib,ib1) 1576 !!! print'(A,5I4,3(2E15.5,3x))', "P1",myproc,is,ikpt,ib,ib1,ff(:) 1577 ! endif 1578 ! !endif 1579 ! enddo 1580 !! call int_fct(ff,ib==ib1,nspinor,paw_dmft,integral) ! here, option==1 even if nspinor==2 1581 !! green%occup%ks(is,ikpt,ib,ib1)=integral 1582 !!! print'(A,5I4,3(2E15.5,3x))', "PP",myproc,is,ikpt,ib,ib1,ff(:) 1583 ! call int_fct(ff,ib==ib1,2,paw_dmft,integral,green%procb(:,ikpt),myproc) ! here, option==1 even if nspinor==2 1584 !!020212 call int_fct(ff,ib==ib1,2,paw_dmft,integral) ! here, option==1 even if nspinor==2 1585 ! green%occup%ks(is,ikpt,ib,ib1)=integral 1586 !!! print'(A,5I4,2E15.8)', "P2",myproc,is,ikpt,ib,ib1,integral 1587 ! ! write(std_out,*) "integral",ikpt,green%occup%ks(is,ikpt,ib,ib1) 1588 ! ! endif 1589 !! write(std_out,'(a,4i6,e14.5,e14.5)') "ks",is,ikpt,ib,ib1,integral 1590 ! enddo ! ib1 1591 ! enddo ! ib 1592 !!020212 endif 1593 ! enddo ! ikpt 1594 ! enddo ! isppol 1595 1596 1597 ! call xmpi_barrier(spacecomm) 1598 ! call xmpi_sum(green%occup%ks,spacecomm,ierr) 1599 1600 1601 !! do is = 1 , nsppol 1602 !! do ikpt = 1, nkpt 1603 !! do ib = 1, mbandc 1604 !! do ib1 = 1, mbandc 1605 !! write(6,'(A,5I4,2E15.8)') "AAAFTERXSUM",myproc,is,ikpt,ib,ib1,green%occup%ks(is,ikpt,ib,ib1) 1606 !! print '(A,5I4,2E15.8)', "AFTERXSUM P",myproc,is,ikpt,ib,ib1,green%occup%ks(is,ikpt,ib,ib1) 1607 !! enddo ! ib1 1608 !! enddo ! ib 1609 !! enddo ! ikpt 1610 !! enddo ! isppol 1611 occup_temp%ks=green%occup%ks 1612 !write(std_out,*) "occup%ks ik1",green%occup%ks(1,1,1,3) 1613 !write(std_out,*) "occup%ks ik2",green%occup%ks(1,2,1,3) 1614 ! - Post-treatment for summation over negative and positive frequencies: 1615 ! necessary in the case of nspinor==2, but valid everywhere 1616 ! N(k,n_1,n_2)= (N(k,n_1,n_2)+ N*(k,n_2,n_1))/2 1617 ! because [G_{k}^{n_1,n_2}(iomega_n)]*= G_{k}^{n_2,n_1}(-iomega_n) 1618 do ib1 = 1, mbandc 1619 do ib = 1, mbandc 1620 do ikpt = 1, nkpt 1621 do is = 1 , nsppol 1622 green%occup%ks(is,ikpt,ib,ib1)=& 1623 & ( occup_temp%ks(is,ikpt,ib,ib1)+ & 1624 & conjg(occup_temp%ks(is,ikpt,ib1,ib)) )/two 1625 enddo ! ib1 1626 enddo ! ib 1627 enddo ! ikpt 1628 enddo ! isppol 1629 !write(std_out,*) "occup%ks ik1 BB",green%occup%ks(1,1,1,3) 1630 !write(std_out,*) "occup%ks ik2 AA",green%occup%ks(1,2,1,3) 1631 call destroy_oper(occup_temp) 1632 do ib1 = 1, mbandc 1633 do ib = 1, mbandc 1634 do ikpt = 1, nkpt 1635 do is = 1 , nsppol 1636 paw_dmft%occnd(1,paw_dmft%include_bands(ib),& 1637 & paw_dmft%include_bands(ib1),ikpt,is)=dreal(green%occup%ks(is,ikpt,ib,ib1)) 1638 paw_dmft%occnd(2,paw_dmft%include_bands(ib),& 1639 & paw_dmft%include_bands(ib1),ikpt,is)=dimag(green%occup%ks(is,ikpt,ib,ib1)) 1640 if(nspinor==1 .and. nsppol==1) then 1641 paw_dmft%occnd(1,paw_dmft%include_bands(ib),& 1642 & paw_dmft%include_bands(ib1),ikpt,is)=two*dreal(green%occup%ks(is,ikpt,ib,ib1)) 1643 paw_dmft%occnd(2,paw_dmft%include_bands(ib),& 1644 & paw_dmft%include_bands(ib1),ikpt,is)=two*dimag(green%occup%ks(is,ikpt,ib,ib1)) 1645 endif 1646 enddo 1647 enddo 1648 enddo 1649 enddo 1650 1651 if(optksloc>0) then 1652 ! - Compute local occupations 1653 call loc_oper(green%occup,paw_dmft,1) 1654 if(abs(prtopt)>2) then 1655 write(message,'(a,a,i10,a)') ch10,& 1656 & " = green%occup%matlu from projection of int(gks(w)) without symetrization" 1657 call wrtout(std_out,message,'COLL') 1658 call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=-1) 1659 endif 1660 1661 ! - Symetrise: continue sum over k-point: Full BZ 1662 call sym_matlu(cryst_struc,green%occup%matlu,pawang,paw_dmft) 1663 if(abs(prtopt)>=2) then 1664 ! write(message,'(a,a,i10,a)') ch10,& 1665 !& " = green%occup%matlu from projection of int(gks(w)) with symetrization" 1666 write(message,'(a,a,i10,a)') ch10,& 1667 & " = Occupation matrix from KS occupations" 1668 call wrtout(std_out,message,'COLL') 1669 call print_matlu(green%occup%matlu,natom,prtopt=3,opt_diag=0) 1670 endif 1671 1672 ! - If the trace of occup matrix in the LOCAL basis was not done 1673 ! before because lpsichiortho/=1 , do it now 1674 if(paw_dmft%lpsichiortho/=1) then 1675 if(abs(prtopt)>0) then 1676 call trace_oper(green%occup,green%charge_ks,green%charge_matlu,2) 1677 green%has_charge_matlu=2 1678 icomp_chloc=1 1679 endif 1680 endif 1681 endif ! optksloc>0: only diagonal elements of G\nunu' are computed 1682 1683 ! - Compute trace over ks density matrix 1684 call trace_oper(green%occup,green%charge_ks,green%charge_matlu,1) 1685 if(abs(prtopt)>0) then 1686 write(message,'(a,a,f12.6)') ch10,& 1687 & " == Total number of electrons from KS green function is :", green%charge_ks 1688 call wrtout(std_out,message,'COLL') 1689 write(message,'(8x,a,f12.6,a)') " (should be",paw_dmft%nelectval,")" 1690 call wrtout(std_out,message,'COLL') 1691 endif 1692 endif ! optksloc 1693 ! ================================================= 1694 1695 1696 ! ================================================= 1697 ! Tests and compute precision on local charge 1698 ! ================================================= 1699 ! - Check that if, renormalized psichi are used, occupations matrices 1700 ! obtained directly from local green function or, through kohn sham 1701 ! occupations are the same. 1702 if(abs(optksloc)==3) then ! optksloc= 3 1703 if(paw_dmft%lpsichiortho==1) then 1704 call diff_matlu("Local_projection_of_kohnsham_occupations ",& 1705 & "Integration_of_local_green_function ",& 1706 & green%occup%matlu,matlu_temp,natom,option,tol4) 1707 write(message,'(2a)') ch10,& 1708 & ' ***** => Calculations of Green function in KS and local spaces are coherent ****' 1709 call wrtout(std_out,message,'COLL') 1710 endif 1711 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-2024 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 DFT. dft_occup mpi_enreg=information about MPI parallelization paw_dmft = data for self-consistent DFT+DMFT calculations.
OUTPUT
paw_dmft = data for self-consistent DFT+DMFT calculations.
SOURCE
3476 subroutine local_ks_green(green,paw_dmft,prtopt) 3477 3478 !Arguments ------------------------------------ 3479 !scalars 3480 type(green_type), intent(in) :: green 3481 type(paw_dmft_type), intent(in) :: paw_dmft 3482 integer, intent(in) :: prtopt 3483 3484 !Local variables ------------------------------ 3485 character(len=500) :: message 3486 integer :: iband,ifreq,ikpt,isppol,itau,lsub,ltau,mbandc,nkpt,nsppol 3487 character(len=1) :: tag_is 3488 character(len=fnlen) :: tmpfil 3489 integer,allocatable :: unitgreenlocks_arr(:) 3490 real(dp) :: beta 3491 real(dp), allocatable :: tau(:) 3492 complex(dpc), allocatable :: loc_ks(:,:,:) 3493 complex(dpc), allocatable :: loc_ks_tau(:,:,:),fw(:),ft(:) 3494 !scalars 3495 !************************************************************************ 3496 mbandc=paw_dmft%mbandc 3497 nkpt=paw_dmft%nkpt 3498 nsppol=paw_dmft%nsppol 3499 ltau=128 3500 ABI_MALLOC(tau,(ltau)) 3501 do itau=1,ltau 3502 tau(itau)=float(itau-1)/float(ltau)/paw_dmft%temp 3503 end do 3504 beta=one/paw_dmft%temp 3505 3506 !Only imaginary frequencies here 3507 if(green%w_type=="real") then 3508 message = ' compute_energy not implemented for real frequency' 3509 ABI_BUG(message) 3510 end if 3511 3512 !========================================= 3513 !Compute local band ks green function 3514 ! should be computed in compute_green: it would be less costly in memory. 3515 !========================================= 3516 ABI_MALLOC(loc_ks,(nsppol,mbandc,paw_dmft%dmft_nwlo)) 3517 if(green%oper(1)%has_operks==1) then 3518 loc_ks(:,:,:)=czero 3519 do isppol=1,nsppol 3520 do iband=1,mbandc 3521 do ifreq=1,paw_dmft%dmft_nwlo 3522 do ikpt=1,nkpt 3523 loc_ks(isppol,iband,ifreq)=loc_ks(isppol,iband,ifreq)+ & 3524 & green%oper(ifreq)%ks(isppol,ikpt,iband,iband)*paw_dmft%wtk(ikpt) 3525 end do 3526 end do 3527 end do 3528 end do 3529 else 3530 message = ' green fct is not computed in ks space' 3531 ABI_BUG(message) 3532 end if 3533 3534 !========================================= 3535 !Compute fourier transformation 3536 !========================================= 3537 3538 ABI_MALLOC(loc_ks_tau,(nsppol,mbandc,ltau)) 3539 ABI_MALLOC(fw,(paw_dmft%dmft_nwlo)) 3540 ABI_MALLOC(ft,(ltau)) 3541 loc_ks_tau(:,:,:)=czero 3542 do isppol=1,nsppol 3543 do iband=1,mbandc 3544 do ifreq=1,paw_dmft%dmft_nwlo 3545 fw(ifreq)=loc_ks(isppol,iband,ifreq) 3546 end do 3547 call fourier_fct(fw,ft,.true.,ltau,-1,paw_dmft) ! inverse fourier 3548 do itau=1,ltau 3549 loc_ks_tau(isppol,iband,itau)=ft(itau) 3550 end do 3551 end do 3552 end do 3553 ABI_FREE(fw) 3554 ABI_FREE(ft) 3555 do isppol=1,nsppol 3556 do iband=1,mbandc 3557 do itau=1,ltau 3558 loc_ks_tau(isppol,iband,itau)=(loc_ks_tau(isppol,iband,itau)+conjg(loc_ks_tau(isppol,iband,itau)))/two 3559 end do 3560 end do 3561 end do 3562 3563 !========================================= 3564 !Print out ksloc green function 3565 !========================================= 3566 if(abs(prtopt)==1) then 3567 ABI_MALLOC(unitgreenlocks_arr,(nsppol)) 3568 do isppol=1,nsppol 3569 write(tag_is,'(i1)')isppol 3570 tmpfil = trim(paw_dmft%filapp)//'Gtau_locks_isppol'//tag_is 3571 write(message,'(3a)') ch10," == Print green function on file ",tmpfil 3572 call wrtout(std_out,message,'COLL') 3573 unitgreenlocks_arr(isppol)=500+isppol-1 3574 open (unit=unitgreenlocks_arr(isppol),file=trim(tmpfil),status='unknown',form='formatted') 3575 rewind(unitgreenlocks_arr(isppol)) 3576 write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenlocks_arr(isppol) 3577 write(message,'(a,a)') ch10,"# New record : First 40 bands" 3578 call wrtout(unitgreenlocks_arr(isppol),message,'COLL') 3579 do lsub=1,mbandc/40+1 3580 do itau=1, ltau 3581 write(message,'(2x,50(e10.3,2x))') tau(itau), & 3582 & (real(loc_ks_tau(isppol,iband,itau)),iband=40*(lsub-1)+1,min(40*lsub,mbandc)) 3583 call wrtout(unitgreenlocks_arr(isppol),message,'COLL') 3584 end do 3585 write(message,'(2x,50(e10.3,2x))') beta, & 3586 & ((-one-real(loc_ks_tau(isppol,iband,1))),iband=40*(lsub-1)+1,min(40*lsub,mbandc)) 3587 call wrtout(unitgreenlocks_arr(isppol),message,'COLL') 3588 if(40*lsub<mbandc) then 3589 write(message,'(a,a,i5,a,i5)') & 3590 & ch10,"# Same record, Following bands : From ", & 3591 & 40*(lsub)," to ",min(40*(lsub+1),mbandc) 3592 call wrtout(unitgreenlocks_arr(isppol),message,'COLL') 3593 end if 3594 end do 3595 ! call flush(unitgreenlocks_arr(isppol)) 3596 end do 3597 ABI_FREE(unitgreenlocks_arr) 3598 end if 3599 3600 !Deallocations 3601 ABI_FREE(loc_ks) 3602 ABI_FREE(loc_ks_tau) 3603 ABI_FREE(tau) 3604 3605 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-2024 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
SOURCE
3125 subroutine newton(cryst_struc,green,paw_dmft,pawang,self,& 3126 & x_input,x_precision,max_iter,f_precision,ierr_hh,opt_noninter,opt_algo) 3127 3128 !Arguments ------------------------------------ 3129 !scalars 3130 type(crystal_t),intent(in) :: cryst_struc 3131 type(green_type),intent(inout) :: green 3132 !type(MPI_type), intent(in) :: mpi_enreg 3133 type(paw_dmft_type), intent(inout) :: paw_dmft 3134 type(pawang_type),intent(in) :: pawang 3135 type(self_type), intent(inout) :: self 3136 integer,intent(in) :: opt_noninter,max_iter 3137 integer,intent(out) :: ierr_hh 3138 real(dp),intent(inout) :: x_input,x_precision 3139 real(dp),intent(inout) :: f_precision 3140 real(dp),intent(in), optional :: opt_algo 3141 !Local variables------------------------------- 3142 integer iter 3143 real(dp) Fx,Fxprime,Fxdouble,xold,x_before,Fxoptimum 3144 real(dp) :: nb_elec_x 3145 integer option,optalgo 3146 logical l_minus,l_plus 3147 real(dp) :: x_minus,x_plus,x_optimum 3148 character(len=500) :: message 3149 ! ********************************************************************* 3150 x_minus=-10_dp 3151 x_plus=-11_dp 3152 xold=-12_dp 3153 3154 if(present(opt_algo)) then 3155 optalgo=opt_algo 3156 else 3157 optalgo=1 ! newton 3158 end if 3159 3160 x_input=paw_dmft%fermie 3161 ierr_hh=0 3162 option =2 ! Halley method 3163 option =1 ! Newton method 3164 3165 !write(std_out,*) "ldaprint",opt_noninter 3166 3167 !--- Start of iterations 3168 write(message,'(a,3a)') " Fermi level Charge Difference" 3169 call wrtout(std_out,message,'COLL') 3170 !do iter=1, 40 3171 !x_input=float(iter)/100_dp 3172 !call function_and_deriv(cryst_struc,f_precision,green,iter,mpi_enreg,paw_dmft,pawang,self,& 3173 !& x_input,x_before,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option) 3174 !write(std_out,*) x_input,Fx 3175 !enddo 3176 !call abi_abort('COLL') 3177 3178 l_minus=.false. 3179 l_plus=.false. 3180 Fxoptimum=1_dp 3181 !======================================== 3182 !start iteration to find fermi level 3183 !======================================== 3184 do iter=1, max_iter 3185 3186 ! ======================================== 3187 ! If zero is located between two values: apply newton method or dichotomy 3188 ! ======================================== 3189 if(l_minus.and.l_plus) then 3190 3191 ! ============================================== 3192 ! Compute the function and derivatives for newton 3193 ! ============================================== 3194 call function_and_deriv(cryst_struc,f_precision,green,iter,paw_dmft,pawang,self,& 3195 & x_input,x_before,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option) 3196 3197 ! Apply stop criterion on Fx 3198 if(abs(Fx) < f_precision) then 3199 ! write(message,'(a,2f12.6)') "Fx,f_precision",Fx,f_precision 3200 ! call wrtout(std_out,message,'COLL') 3201 x_precision=x_input-xold 3202 return 3203 end if 3204 if(iter==max_iter) then 3205 write(message,'(a,2f12.6)') " Fermi level could not be found" 3206 call wrtout(std_out,message,'COLL') 3207 x_input=x_optimum 3208 ierr_hh=-123 3209 return 3210 end if 3211 3212 ! Cannot divide by Fxprime if too small 3213 if(abs(Fxprime) .le. 1.e-15)then 3214 ierr_hh=-314 3215 write(message,'(a,f12.7)') "Fxprime=",Fxprime 3216 call wrtout(std_out,message,'COLL') 3217 return 3218 end if 3219 3220 x_precision=x_input-xold 3221 3222 ! ============================================== 3223 ! Newton/Halley's formula for next iteration 3224 ! ============================================== 3225 xold=x_input 3226 if(option==1) x_input=x_input-Fx/Fxprime 3227 if(option==2) x_input=x_input-2*Fx*Fxprime/(2*Fxprime**2-Fx*Fxdouble) 3228 3229 ! ============================================== 3230 ! If newton does not work well, use dichotomy. 3231 ! ============================================== 3232 if(x_input<x_minus.or.x_input>x_plus) then 3233 call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,& 3234 & Fx,opt_noninter,nb_elec_x,xold) 3235 write(message,'(a,3f12.6)') " ---",x_input,nb_elec_x,Fx 3236 call wrtout(std_out,message,'COLL') 3237 if(Fx>0) then 3238 x_plus=xold 3239 else if(Fx<0) then 3240 x_minus=xold 3241 end if 3242 x_input=(x_plus+x_minus)/2.d0 3243 3244 end if 3245 ! write(std_out,'(a,2f12.6)') " Q(xold) and dQ/dx=",Fx,Fxprime 3246 ! write(std_out,'(a,f12.6)') " => new Fermi level",x_input 3247 ! ======================================== 3248 ! Locate zero between two values 3249 ! ======================================== 3250 else 3251 call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,& 3252 & Fx,opt_noninter,nb_elec_x,x_input) 3253 write(message,'(a,3f12.6)') " --",x_input,nb_elec_x,Fx 3254 ! Possible improvement for large systems, removed temporarely for 3255 ! automatic tests: more study is necessary: might worsen the convergency 3256 ! if(iter==1) then 3257 ! f_precision=max(abs(Fx/50),f_precision) 3258 ! write(message,'(a,4x,a,e12.6)') ch10," Precision required changed to:",f_precision 3259 ! call wrtout(std_out,message,'COLL') 3260 ! endif 3261 call wrtout(std_out,message,'COLL') 3262 if(Fx>0) then 3263 l_plus=.true. 3264 x_plus=x_input 3265 x_input=x_input-0.02 3266 else if(Fx<0) then 3267 l_minus=.true. 3268 x_minus=x_input 3269 x_input=x_input+0.02 3270 end if 3271 3272 end if 3273 3274 if(abs(Fx)<abs(Fxoptimum)) then 3275 Fxoptimum=Fx 3276 x_optimum=x_input 3277 end if 3278 3279 3280 3281 ! if(myid==master) then 3282 ! write(std_out,'(a,i4,3f12.6)') "i,xnew,F,Fprime",i,x_input,Fx,Fxprime 3283 ! endif 3284 3285 3286 ! ! Apply stop criterion on x 3287 ! if(abs(x_input-xold) .le. x_input*x_precision) then 3288 ! ! write(std_out,'(a,4f12.6)') "x_input-xold, x_precision*x_input "& 3289 ! !& ,x_input-xold,x_precision*x_input,x_precision 3290 ! f_precision=Fx 3291 ! return 3292 ! endif 3293 3294 end do 3295 !--- End of iterations 3296 3297 3298 ierr_hh=1 3299 return 3300 3301 CONTAINS !======================================================================================== 3302 !-----------------------------------------------------------------------
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
SOURCE
2698 subroutine occup_green_tau(green) 2699 2700 !Arguments ------------------------------------ 2701 !type 2702 type(green_type),intent(inout) :: green 2703 2704 !local variables------------------------------- 2705 integer :: iatom,lpawu 2706 complex(dpc), allocatable :: shift(:) 2707 ! ********************************************************************* 2708 2709 ABI_MALLOC(shift,(green%oper_tau(1)%natom)) 2710 do iatom=1,green%oper_tau(1)%natom 2711 shift(iatom)=-cone 2712 lpawu=green%oper_tau(1)%matlu(iatom)%lpawu 2713 if(lpawu/=-1) then 2714 ! do isppol=1,green%oper_tau(1)%nsppol 2715 green%occup_tau%matlu(iatom)%mat(:,:,:,:,:)= & 2716 & green%oper_tau(1)%matlu(iatom)%mat(:,:,:,:,:) 2717 endif 2718 enddo 2719 call shift_matlu(green%occup_tau%matlu,green%occup_tau%natom,shift) 2720 ABI_FREE(shift) 2721 2722 2723 end subroutine occup_green_tau
m_green/occupfd [ Functions ]
[ Top ] [ m_green ] [ Functions ]
NAME
occupfd
FUNCTION
INPUTS
OUTPUT
SOURCE
2739 function occupfd(eig,fermie,temp) 2740 2741 2742 !Arguments ------------------------------------ 2743 !type 2744 ! Integrate analytic tail 1/(iw-mu) 2745 real(dp),intent(in) :: eig,fermie,temp 2746 real(dp) :: occupfd 2747 !Local variables------------------------------- 2748 ! ********************************************************************* 2749 2750 if((eig-fermie) > zero) then 2751 occupfd=exp(-(eig-fermie)/temp)/(one+exp(-(eig-fermie)/temp)) 2752 else 2753 occupfd=one/(one+exp((eig-fermie)/temp)) 2754 endif 2755 2756 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" 5 print k-resolved 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
SOURCE
635 subroutine print_green(char1,green,option,paw_dmft,pawprtvol,opt_wt,opt_decim) 636 637 !Arguments ------------------------------------ 638 !type 639 type(paw_dmft_type), intent(in) :: paw_dmft 640 type(green_type),intent(inout) :: green 641 integer,intent(in) :: option,pawprtvol 642 integer,intent(in),optional :: opt_wt,opt_decim 643 character(len=*), intent(in) :: char1 644 645 !local variables------------------------------- 646 integer :: iall,iatom,ib,ifreq,ikpt,im,ispinor,isppol,itau 647 integer :: lsub,mbandc,natom,ndim,nkpt,nspinor,nsppol,optwt,spf_unt,spfkresolved_unt,spcorb_unt 648 character(len=2000) :: message 649 integer,allocatable :: unitgreenfunc_arr(:) 650 integer,allocatable :: unitgreenloc_arr(:) 651 character(len=fnlen) :: tmpfil 652 character(len=1) :: tag_is,tag_is2 653 character(len=10) :: tag_at 654 character(len=3) :: tag_ik 655 real(dp) :: re, ima 656 complex(dpc), allocatable :: sf(:),sf_corr(:) 657 ! ********************************************************************* 658 if(present(opt_wt)) then 659 optwt=opt_wt 660 else 661 optwt=1 662 endif 663 664 665 if(pawprtvol>200) then 666 endif 667 natom=green%oper(1)%natom 668 nsppol=paw_dmft%nsppol 669 nspinor=paw_dmft%nspinor 670 mbandc=paw_dmft%mbandc 671 nkpt=paw_dmft%nkpt 672 673 ! == Print local Green Function 674 if(option==1.or.option==3) then 675 ABI_MALLOC(unitgreenfunc_arr,(natom*nsppol*nspinor)) 676 iall=0 677 do iatom=1,natom 678 if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then 679 ndim=2*green%oper(1)%matlu(iatom)%lpawu+1 680 call int2char4(iatom,tag_at) 681 ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!') 682 do isppol=1,nsppol 683 write(tag_is,'(i1)')isppol 684 do ispinor=1,nspinor 685 iall=iall+1 686 write(tag_is2,'(i1)')ispinor 687 688 ! == Create names 689 if(optwt==1) then 690 tmpfil = & 691 & trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-omega_iatom'// & 692 & trim(tag_at)//'_isppol'//tag_is//'_ispinor'//tag_is2 693 else 694 tmpfil = & 695 & trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-tau_iatom'// & 696 & trim(tag_at)//'_isppol'//tag_is//'_ispinor'//tag_is2 697 endif 698 if(iall<=4) then 699 write(message,'(3a)') ch10," == Print green function on file ",tmpfil 700 call wrtout(std_out,message,'COLL') 701 elseif(iall==5) then 702 write(message,'(3a)') ch10," == following values are printed in files" 703 call wrtout(std_out,message,'COLL') 704 endif 705 unitgreenfunc_arr(iall)=300+iall-1 706 if(optwt==1.or.green%fileprt_tau==0) then 707 open(unit=unitgreenfunc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted',position='append') 708 else if(optwt==2.and.green%fileprt_tau==1) then 709 open(unit=unitgreenfunc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted',position='append') 710 endif 711 712 ! == Write in files 713 ! rewind(unitgreenfunc_arr(iall)) 714 ! write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenfunc_arr(iall) 715 ! call wrtout(std_out,message,'COLL') 716 write(message,'(a,a)') ch10,"# New record :" 717 call wrtout(unitgreenfunc_arr(iall),message,'COLL') 718 if(optwt==1) then 719 do ifreq=1,green%nw 720 if(present(opt_decim)) then 721 write(message,'(2x,30(e23.16,2x))') & 722 & green%omega(ifreq),& 723 & (green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor),im=1,ndim) 724 else 725 write(message,'(2x,30(e10.3,2x))') & 726 & green%omega(ifreq),& 727 & (green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor),im=1,ndim) 728 endif 729 call wrtout(unitgreenfunc_arr(iall),message,'COLL') 730 re=real(green%oper(ifreq)%matlu(iatom)%mat(1,1,isppol,ispinor,ispinor)) 731 ima=aimag(green%oper(ifreq)%matlu(iatom)%mat(1,1,isppol,ispinor,ispinor)) 732 ! write(228,*) green%omega(ifreq),re/(re**2+ima**2),ima/(re**2+ima**2)+green%omega(ifreq) 733 enddo 734 else 735 do itau=1,green%dmftqmc_l 736 write(message,'(2x,30(e10.3,2x))') & 737 & green%tau(itau),& 738 & (green%oper_tau(itau)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor),im=1,ndim) 739 call wrtout(unitgreenfunc_arr(iall),message,'COLL') 740 enddo 741 write(message,'(2x,30(e10.3,2x))') & 742 & one/paw_dmft%temp,& 743 & (-green%oper_tau(1)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor)-one,im=1,ndim) 744 call wrtout(unitgreenfunc_arr(iall),message,'COLL') 745 endif 746 enddo 747 enddo ! isppol 748 endif ! lpawu=/-1 749 enddo ! iatom 750 ABI_FREE(unitgreenfunc_arr) 751 endif 752 753 ! == Print ks green function 754 if((option==2.or.option==3).and.green%oper(1)%has_operks==1) then 755 ABI_MALLOC(unitgreenloc_arr,(nsppol*nkpt)) 756 iall=0 757 do isppol = 1 , nsppol 758 write(tag_is,'(i1)')isppol 759 do ikpt = 1, nkpt 760 write(tag_ik,'(i3)')ikpt 761 ! do ib1 = 1, mbandc 762 iall=iall+1 763 764 ! == Create names 765 if(optwt==1) then 766 tmpfil = trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-omega_isppol'//tag_is//'_ikpt'//trim(adjustl(tag_ik)) 767 else 768 tmpfil = trim(paw_dmft%filapp)//'Green-'//trim(char1)//'-tau_isppol'//tag_is//'_ikpt'//trim(adjustl(tag_ik)) 769 endif 770 if(iall<=4) then 771 write(message,'(3a)') ch10," == Print green function on file ",tmpfil 772 call wrtout(std_out,message,'COLL') 773 elseif(iall==5) then 774 write(message,'(3a)') ch10," == following values are printed in files" 775 call wrtout(std_out,message,'COLL') 776 endif 777 unitgreenloc_arr(iall)=400+iall-1 778 open (unit=unitgreenloc_arr(iall),file=trim(tmpfil),status='unknown',form='formatted') 779 ! rewind(unitgreenloc_arr(iall)) 780 ! write(message,'(a,a,a,i4)') 'opened file : ', trim(tmpfil), ' unit', unitgreenloc_arr(iall) 781 ! call wrtout(std_out,message,'COLL') 782 783 ! == Write in files 784 write(message,'(a,a)') ch10,"# New record : First 20 bands" 785 call wrtout(unitgreenloc_arr(iall),message,'COLL') 786 ! call flush(std_out) 787 do lsub=1,mbandc/20+1 788 if(optwt==1) then 789 do ifreq=1,green%nw 790 ! call flush(std_out) 791 write(message,'(2x,50(e10.3,2x))') & 792 & green%omega(ifreq), & 793 & (green%oper(ifreq)%ks(isppol,ikpt,ib,ib),ib=20*(lsub-1)+1,min(20*lsub,mbandc)) 794 call wrtout(unitgreenloc_arr(iall),message,'COLL') 795 re=real(green%oper(ifreq)%ks(isppol,ikpt,1,1)) 796 ima=aimag(green%oper(ifreq)%ks(isppol,ikpt,1,1)) 797 if(ikpt==1) write(229,*) green%omega(ifreq),re/(re**2+ima**2),ima/(re**2+ima**2)+green%omega(ifreq) 798 enddo 799 else 800 if(green%use_oper_tau_ks==1) then 801 do itau=1,green%dmftqmc_l 802 ! call flush(std_out) 803 write(message,'(2x,50(e10.3,2x))') & 804 & green%tau(itau), & 805 & (green%oper_tau(itau)%ks(isppol,ikpt,ib,ib),ib=20*(lsub-1)+1,min(20*lsub,mbandc)) 806 call wrtout(unitgreenloc_arr(iall),message,'COLL') 807 enddo 808 endif 809 endif 810 if(20*lsub<mbandc) write(message,'(a,a,i5,a,i5)') & 811 & ch10,"# Same record, Following bands : From ", & 812 & 20*(lsub)," to ",min(20*(lsub+1),mbandc) 813 call wrtout(unitgreenloc_arr(iall),message,'COLL') 814 enddo 815 enddo ! ikpt 816 enddo ! isppol 817 ABI_FREE(unitgreenloc_arr) 818 endif 819 820 if((green%w_type=="real".and.option>=4).and.green%oper(1)%has_operks==1) then 821 write(message,'(a,a)') ch10," == About to print spectral function" 822 call wrtout(std_out,message,'COLL') 823 if (option==4) then 824 tmpfil = trim(paw_dmft%filapp)//'SpFunc-'//trim(char1) 825 if (open_file(tmpfil, message, newunit=spf_unt, status='unknown', form='formatted') /= 0) then 826 ABI_ERROR(message) 827 end if 828 endif 829 if (option==5) then 830 tmpfil = trim(paw_dmft%filapp)//'_DFTDMFT_SpectralFunction_kresolved_'//trim(char1) 831 if (open_file(tmpfil, message, newunit=spfkresolved_unt, status='unknown', form='formatted') /= 0) then 832 ABI_ERROR(message) 833 end if 834 endif 835 ABI_MALLOC(sf,(green%nw)) 836 ABI_MALLOC(sf_corr,(green%nw)) 837 iall=0 838 if (option==5) then 839 do ikpt = 1, nkpt 840 sf=czero 841 do ifreq=1,green%nw 842 do isppol = 1 , nsppol 843 do ib=1,mbandc 844 sf(ifreq)=sf(ifreq)+green%oper(ifreq)%ks(isppol,ikpt,ib,ib) 845 enddo 846 enddo 847 write(message,*) green%omega(ifreq)*Ha_eV,(-aimag(sf(ifreq)))/pi/Ha_eV,ikpt 848 call wrtout(spfkresolved_unt,message,'COLL') 849 enddo 850 write(message,*) 851 call wrtout(spfkresolved_unt,message,'COLL') 852 enddo 853 write(message,*) ch10 854 call wrtout(spfkresolved_unt,message,'COLL') 855 ! 856 ! do isppol = 1 , nsppol 857 ! do ikpt = 1, nkpt 858 ! do ib=1,mbandc 859 ! sf=czero 860 ! write(71,*) 861 ! write(71,*) "#", ikpt, ib 862 ! do ifreq=1,green%nw 863 ! sf(ifreq)=sf(ifreq)+green%oper(ifreq)%ks(isppol,ikpt,ib,ib) 864 ! write(71,*) green%omega(ifreq)*Ha_eV,(-aimag(sf(ifreq)))/pi/Ha_eV,ikpt 865 ! enddo 866 ! enddo 867 ! enddo 868 ! enddo 869 endif 870 if (option==4) then 871 sf=czero 872 do isppol = 1 , nsppol 873 do ikpt = 1, nkpt 874 do ib=1,mbandc 875 do ifreq=1,green%nw 876 sf(ifreq)=sf(ifreq)+green%oper(ifreq)%ks(isppol,ikpt,ib,ib)*green%oper(1)%wtk(ikpt) 877 enddo 878 enddo 879 enddo 880 enddo 881 endif 882 if(paw_dmft%dmft_kspectralfunc==1) then 883 do iatom=1,natom 884 if(green%oper(1)%matlu(iatom)%lpawu.ne.-1) then 885 sf_corr=czero 886 call int2char4(iatom,tag_at) 887 ABI_CHECK((tag_at(1:1)/='#'),'Bug: string length too short!') 888 tmpfil = trim(paw_dmft%filapp)//'_DFTDMFT_spectralfunction_orb_'//trim(char1)//'_iatom'//trim(tag_at) 889 if (open_file(tmpfil, message, newunit=spcorb_unt, status='unknown', form='formatted') /= 0) then 890 ABI_ERROR(message) 891 end if 892 write(message,*) "#", nspinor,nsppol,ndim,green%nw 893 call wrtout(spcorb_unt,message,'COLL') 894 write(message,*) "#", green%oper(1)%matlu(iatom)%lpawu 895 call wrtout(spcorb_unt,message,'COLL') 896 ndim=2*green%oper(1)%matlu(iatom)%lpawu+1 897 do isppol = 1 , nsppol 898 do ispinor=1, nspinor 899 do im=1,ndim 900 do ifreq=1,green%nw 901 sf_corr(ifreq)=sf_corr(ifreq)+ green%oper(ifreq)%matlu(iatom)%mat(im,im,isppol,ispinor,ispinor) 902 enddo 903 enddo 904 enddo 905 enddo 906 do ifreq=1,green%nw 907 write(message,*) green%omega(ifreq),(-aimag(sf_corr(ifreq)))/3.141592653589793238_dp 908 call wrtout(spcorb_unt,message,'COLL') 909 enddo 910 close(spcorb_unt) 911 endif 912 enddo 913 endif 914 if (option==4) then 915 do ifreq=1,green%nw 916 write(message,*) green%omega(ifreq),(-aimag(sf(ifreq)))/3.141592653589793238_dp 917 call wrtout(spf_unt,message,'COLL') 918 enddo 919 endif 920 close(spf_unt) 921 ABI_FREE(sf) 922 ABI_FREE(sf_corr) 923 endif 924 925 if(optwt==2.and.(option==1.or.option==3)) then 926 green%fileprt_tau=1 ! file for G(tau) has been created here 927 endif 928 929 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
SOURCE
535 subroutine printocc_green(green,option,paw_dmft,pawprtvol,opt_weissgreen,chtype) 536 537 !Arguments ------------------------------------ 538 !type 539 type(paw_dmft_type), intent(in) :: paw_dmft 540 type(green_type),intent(in) :: green 541 integer,intent(in) :: option,pawprtvol 542 integer,optional,intent(in) :: opt_weissgreen 543 character(len=*), optional, intent(in) :: chtype 544 545 !local variables------------------------------- 546 character(len=500) :: message 547 integer :: optweissgreen,i_tau 548 ! ********************************************************************* 549 if(present(opt_weissgreen)) then 550 optweissgreen=opt_weissgreen 551 else 552 optweissgreen=2 553 endif 554 555 556 if(mod(option,4)==1) then 557 if(optweissgreen==2) then 558 if(present(chtype)) then 559 write(message,'(4a)') ch10," == The ",trim(chtype)," occupations are == " 560 else 561 write(message,'(2a)') ch10," == The occupations (integral of the Green function) are == " 562 endif 563 else if(optweissgreen==1) then 564 write(message,'(2a)') ch10," == The integrals of the Weiss function are == " 565 endif 566 call wrtout(std_out,message,'COLL') 567 call print_oper(green%occup,option,paw_dmft,pawprtvol) 568 endif 569 if(mod(option,4)>=2) then 570 if(optweissgreen==2) then 571 write(message,'(2a)') ch10," == The occupations (value of G(tau) for tau=0-) are == " 572 else if(optweissgreen==1) then 573 write(message,'(2a)') ch10," == Values of G_0(tau) for tau=0- are == " 574 endif 575 call wrtout(std_out,message,'COLL') 576 call print_oper(green%occup_tau,option,paw_dmft,pawprtvol) 577 ! write(message,'(2a)') ch10," == check: occupations from Green functions are == " 578 ! call wrtout(std_out,message,'COLL') 579 ! call print_oper(green%occup,1,paw_dmft,pawprtvol) 580 if(mod(option,4)>=3) then 581 call diff_matlu("Local occup from integral of G(w) ",& 582 & "Local occup from G(tau=0-) ",& 583 & green%occup%matlu,green%occup_tau%matlu,paw_dmft%natom,1,tol4) 584 write(message,'(2a)') ch10,& 585 & ' ***** => Calculations of occupations in omega and tau spaces are coherent ****' 586 call wrtout(std_out,message,'COLL') 587 endif 588 endif 589 590 if(present(chtype)) then 591 if(paw_dmft%prtvol>=4.and.& 592 & (chtype=="DMFT (end of DMFT loop)".or.chtype=="converged DMFT")& 593 & .and.green%occup%has_opermatlu==1) then 594 write(message,'(4a)') ch10," == The DFT+DMFT occupation matrix for correlated electrons is == " 595 call wrtout(ab_out,message,'COLL') 596 call print_matlu(green%occup%matlu,paw_dmft%natom,pawprtvol,opt_ab_out=1) 597 write(message,'(a)') " " 598 call wrtout(ab_out,message,'COLL') 599 endif 600 endif 601 602 if(mod(option,4)>=2) then 603 i_tau=1 604 if(optweissgreen==1) i_tau=-1 605 call trace_matlu(green%occup_tau%matlu,paw_dmft%natom,itau=i_tau) 606 endif 607 608 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-2024 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
SOURCE
2634 subroutine spline_fct(fw1,fw2,opt_spline,paw_dmft) 2635 2636 !Arguments ------------------------------------ 2637 !type 2638 integer,intent(in) :: opt_spline 2639 type(paw_dmft_type), intent(in) :: paw_dmft 2640 complex(dpc), intent(inout) :: fw1(:) 2641 complex(dpc), intent(inout) :: fw2(:) 2642 integer :: size_fw1 2643 integer :: size_fw2 2644 real(dp), allocatable :: omega_li(:) 2645 2646 ! ********************************************************************* 2647 2648 size_fw1 = size(fw1) 2649 size_fw2 = size(fw2) 2650 2651 ! == inverse fourier transform 2652 if(opt_spline==-1) then 2653 2654 ! allocate(fw1(0:paw_dmft%dmft_nwlo-1)) 2655 if(paw_dmft%dmft_log_freq==1) then 2656 ABI_MALLOC(omega_li,(1:size_fw2)) 2657 call construct_nwli_dmft(paw_dmft,size_fw2,omega_li) 2658 call spline_c(size_fw1,size_fw2,paw_dmft%omega_lo(1:size_fw1),& 2659 & omega_li(1:size_fw2),fw2,fw1) 2660 ABI_FREE(omega_li) 2661 else 2662 fw2=fw1 2663 endif 2664 2665 ! == direct fourier transform 2666 else if(opt_spline==1) then 2667 2668 2669 if(paw_dmft%dmft_log_freq==1) then 2670 ABI_MALLOC(omega_li,(1:size_fw2)) 2671 call construct_nwli_dmft(paw_dmft,size_fw2,omega_li) 2672 call spline_c(size_fw2,size_fw1,omega_li(1:size_fw2),& 2673 & paw_dmft%omega_lo(1:size_fw1),fw1,fw2) 2674 ABI_FREE(omega_li) 2675 else 2676 fw1=fw2 2677 endif 2678 2679 endif 2680 2681 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
SOURCE
3415 subroutine compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,& 3416 & Fx,opt_noninter,nb_elec_x,fermie) 3417 3418 3419 3420 !Arguments ------------------------------------ 3421 !scalars 3422 type(crystal_t),intent(in) :: cryst_struc 3423 type(green_type),intent(inout) :: green 3424 !type(MPI_type), intent(in) :: mpi_enreg 3425 type(paw_dmft_type), intent(inout) :: paw_dmft 3426 type(pawang_type),intent(in) :: pawang 3427 type(self_type), intent(inout) :: self 3428 integer,intent(in) :: opt_noninter 3429 real(dp),intent(in) :: fermie 3430 real(dp),intent(out) :: Fx,nb_elec_x 3431 ! ********************************************************************* 3432 paw_dmft%fermie=fermie 3433 call compute_green(cryst_struc,green,paw_dmft,pawang,0,self,opt_self=1,& 3434 & opt_nonxsum=1,opt_nonxsum2=1) 3435 call integrate_green(cryst_struc,green,paw_dmft,pawang,prtopt=0,& 3436 & opt_ksloc=-1) !,opt_nonxsum=1) 3437 ! opt_ksloc=-1, compute total charge 3438 nb_elec_x=green%charge_ks 3439 Fx=nb_elec_x-paw_dmft%nelectval 3440 3441 if(opt_noninter==1) then 3442 end if 3443 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)
SOURCE
3325 subroutine function_and_deriv(cryst_struc,f_precision,green,iter,paw_dmft,pawang,self& 3326 & ,x_input,x_old,x_precision,Fx,Fxprime,Fxdouble,opt_noninter,option) 3327 3328 3329 3330 !Arguments ------------------------------------ 3331 !scalars 3332 type(crystal_t),intent(in) :: cryst_struc 3333 type(green_type),intent(inout) :: green 3334 !type(MPI_type), intent(in) :: mpi_enreg 3335 type(paw_dmft_type), intent(inout) :: paw_dmft 3336 type(pawang_type),intent(in) :: pawang 3337 type(self_type), intent(inout) :: self 3338 integer,intent(in) :: iter,opt_noninter,option 3339 real(dp),intent(inout) :: f_precision,x_input,x_precision 3340 real(dp),intent(out) :: Fx,Fxprime,Fxdouble 3341 real(dp),intent(inout) :: x_old 3342 !Local variables------------------------------- 3343 real(dp) :: deltax,nb_elec_x,Fxmoins,Fxplus,xmoins,xplus,x0 3344 character(len=500) :: message 3345 ! ********************************************************************* 3346 3347 ! choose deltax: for numeric evaluation of derivative 3348 if(iter==1) then 3349 ! deltax=0.02 3350 end if 3351 ! deltax=max((x_input-x_old)/10.d0,min(0.00001_dp,x_precision/100_dp)) 3352 deltax=min(0.00001_dp,x_precision/100_dp) ! small but efficient 3353 ! endif 3354 ! write(std_out,*) "iter,x_input,deltax",iter,x_input,deltax 3355 x0=x_input 3356 xmoins=x0-deltax 3357 xplus=x0+deltax 3358 3359 call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,& 3360 & Fx,opt_noninter,nb_elec_x,x0) 3361 3362 write(message,'(a,3f12.6)') " - ",x0,nb_elec_x,Fx 3363 call wrtout(std_out,message,'COLL') 3364 ! write(std_out,*) "Fx", Fx 3365 if(abs(Fx)<f_precision) return 3366 3367 call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,& 3368 & Fxplus,opt_noninter,nb_elec_x,xplus) 3369 3370 write(message,'(a,3f12.6)') " - ",xplus,nb_elec_x,Fxplus 3371 call wrtout(std_out,message,'COLL') 3372 3373 if(option==2) then 3374 call compute_nb_elec(cryst_struc,green,paw_dmft,pawang,self,& 3375 & Fxmoins,opt_noninter,nb_elec_x,xmoins) 3376 3377 write(message,'(a,3f12.6)') " - ",xmoins,nb_elec_x,Fxmoins 3378 call wrtout(std_out,message,'COLL') 3379 end if 3380 3381 if(option==1) then 3382 Fxprime=(Fxplus-Fx)/deltax 3383 else if (option==2) then 3384 Fxprime=(Fxplus-Fxmoins)/(2*deltax) 3385 Fxdouble=(Fxplus+Fxmoins-2*Fx)/(deltax**2) 3386 end if 3387 ! write(std_out,*) "after computation of Fxprime",myid 3388 if(Fxprime<zero) then 3389 write(message,'(a,f12.6)') " Warning: slope of charge versus fermi level is negative !", Fxprime 3390 call wrtout(std_out,message,'COLL') 3391 end if 3392 x_old=x_input 3393 3394 return 3395 end subroutine function_and_deriv