Line data Source code
1 : MODULE m_kk_cutoff
2 :
3 : USE m_trapz
4 : USE m_types
5 : USE m_juDFT
6 : USE m_constants
7 :
8 : IMPLICIT NONE
9 :
10 : CONTAINS
11 :
12 26 : SUBROUTINE kk_cutoff(im,noco,l_mperp,l,jspins,eMesh,cutoff,scalingFactor)
13 :
14 : !This Subroutine determines the cutoff energy for the kramers-kronig-integration
15 : !This cutoff energy is defined so that the integral over the projDOS up to this cutoff
16 : !is equal to 2*(2l+1) (the number of states in the correlated shell) or not to small
17 :
18 : COMPLEX, INTENT(IN) :: im(:,-lmaxU_const:,-lmaxU_const:,:)
19 : TYPE(t_noco), INTENT(IN) :: noco
20 : LOGICAL, INTENT(IN) :: l_mperp
21 : INTEGER, INTENT(IN) :: l
22 : INTEGER, INTENT(IN) :: jspins
23 : REAL, INTENT(IN) :: eMesh(:)
24 : INTEGER, INTENT(INOUT) :: cutoff(:,:)
25 : REAL, INTENT(INOUT) :: scalingFactor(:)
26 :
27 : INTEGER :: m,ispin,spins_cut,ne
28 : REAL :: lowerBound,upperBound,integral,n_states
29 : REAL :: ec,del,eb,et
30 26 : REAL, ALLOCATABLE :: projDOS(:,:)
31 :
32 26 : ne = SIZE(eMesh)
33 26 : del = eMesh(2)-eMesh(1)
34 26 : eb = eMesh(1)
35 26 : et = eMesh(ne)
36 :
37 :
38 : !Calculate the trace over m,mp of the Imaginary Part matrix to obtain the projected DOS
39 : !n_f(e) = -1/pi * TR[Im(G_f(e))]
40 352956 : ALLOCATE(projDOS(ne,jspins),source=0.0)
41 78 : DO ispin = 1, jspins
42 366 : DO m = -l , l
43 1893940 : projDOS(:,ispin) = projDOS(:,ispin) - 1/pi_const * REAL(im(:,m,m,ispin))
44 : ENDDO
45 : ENDDO
46 :
47 :
48 26 : spins_cut = MERGE(1,jspins,noco%l_noco.AND.l_mperp)
49 26 : n_states = (2*l+1) * MERGE(2.0,2.0/jspins,noco%l_noco.AND.l_mperp)
50 :
51 78 : cutoff(:,1) = 1 !we don't modify the lower bound
52 78 : cutoff(:,2) = ne
53 78 : scalingFactor(:) = 1.0
54 :
55 74 : DO ispin = 1, spins_cut
56 :
57 : !----------------------------------------
58 : !Check the integral up to the hard cutoff
59 : !----------------------------------------
60 10848 : IF(spins_cut.EQ.1 .AND.jspins.EQ.2) projDOS(:,1) = projDOS(:,1) + projDOS(:,2)
61 :
62 : !Initial complete integral
63 48 : integral = trapz(projDOS(:,ispin),del,ne)
64 :
65 : #ifdef CPP_DEBUG
66 : WRITE(*,*) "Integral over DOS: ", integral
67 : #endif
68 :
69 74 : IF(integral.LT.n_states) THEN
70 : !If we are calculating the greens function for a d-band this is expected to happen
71 4 : IF(l.EQ.2) THEN
72 4 : scalingFactor(ispin) = n_states/integral
73 :
74 : #ifdef CPP_DEBUG
75 : WRITE(*,9000) l,ispin,scalingFactor(ispin)
76 : 9000 FORMAT("Scaling the DOS for l=",I1," and spin ",I1," factor: ",f14.8)
77 : #endif
78 4 : IF(scalingFactor(ispin).GT.1.25) CALL juDFT_warn("scaling factor >1.25 -> increase elup(<1htr) or numbands",calledby="kk_cutoff")
79 0 : ELSE IF(integral.LT.n_states-0.1) THEN
80 : ! If the integral is to small we terminate here to avoid problems
81 0 : CALL juDFT_warn("Integral over DOS too small for f -> increase elup(<1htr) or numbands", calledby="kk_cutoff")
82 : ENDIF
83 : ELSE
84 :
85 : !IF the integral is bigger than 2l+1, search for the cutoff using the bisection method
86 : lowerBound = eb
87 : upperBound = et
88 :
89 660 : DO WHILE(ABS(upperBound-lowerBound).GT.del/2.0)
90 :
91 616 : ec = (lowerBound+upperBound)/2.0
92 616 : cutoff(ispin,2) = INT((ec-eb)/del)+1
93 :
94 : !Integrate the DOS up to the cutoff
95 616 : integral = trapz(projDOS(:,ispin),del,cutoff(ispin,2))
96 :
97 660 : IF(ABS(integral-n_states).LT.1e-12) THEN
98 : EXIT
99 616 : ELSE IF(integral.LT.n_states) THEN
100 : !integral to small -> choose the right interval
101 : lowerBound = ec
102 273 : ELSE IF(integral.GT.n_states) THEN
103 : !integral to big -> choose the left interval
104 273 : upperBound = ec
105 : END IF
106 :
107 : ENDDO
108 :
109 : #ifdef CPP_DEBUG
110 : WRITE(*,*) "CALCULATED CUTOFF: ", cutoff(ispin,2)
111 : WRITE(*,*) "CORRESPONDING ENERGY", ec
112 : WRITE(*,*) "INTEGRAL OVER projDOS with cutoff: ", integral
113 : #endif
114 : !Copy cutoff to second spin if only one was calculated
115 44 : IF(spins_cut.EQ.1 .AND. jspins.EQ.2) THEN
116 4 : cutoff(2,2) = cutoff(1,2)
117 4 : scalingFactor(2) = scalingFactor(1)
118 : ENDIF
119 : ENDIF
120 : ENDDO
121 :
122 26 : END SUBROUTINE kk_cutoff
123 :
124 2 : SUBROUTINE kk_cutoffRadial(uu,ud,du,dd,noco,scalarGF,l_mperp,&
125 2 : l,input,eMesh,cutoff,scalingFactor)
126 :
127 : COMPLEX, INTENT(IN) :: uu(:,-lmaxU_const:,-lmaxU_const:,:)
128 : COMPLEX, INTENT(IN) :: ud(:,-lmaxU_const:,-lmaxU_const:,:)
129 : COMPLEX, INTENT(IN) :: du(:,-lmaxU_const:,-lmaxU_const:,:)
130 : COMPLEX, INTENT(IN) :: dd(:,-lmaxU_const:,-lmaxU_const:,:)
131 : TYPE(t_noco), INTENT(IN) :: noco
132 : TYPE(t_scalarGF), INTENT(IN) :: scalarGF
133 : LOGICAL, INTENT(IN) :: l_mperp
134 : INTEGER, INTENT(IN) :: l
135 : TYPE(t_input), INTENT(IN) :: input
136 : REAL, INTENT(IN) :: eMesh(:)
137 : INTEGER, INTENT(INOUT) :: cutoff(:,:)
138 : REAL, INTENT(INOUT) :: scalingFactor(:)
139 :
140 2 : COMPLEX, ALLOCATABLE :: im(:,:,:,:)
141 :
142 : INTEGER :: jspin,m,mp,spin1,spin2
143 :
144 : !calculate the spherical average from the original greens function
145 1058638 : ALLOCATE(im(SIZE(uu,1),-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,SIZE(uu,4)),source=cmplx_0)
146 6 : DO jspin = 1, SIZE(im,4)
147 4 : IF(jspin < 3) THEN
148 4 : spin1 = jspin
149 4 : spin2 = jspin
150 : ELSE
151 : spin1 = 2
152 : spin2 = 1
153 : ENDIF
154 : !$OMP parallel do default(none) &
155 : !$OMP shared(scalarGF,jspin,spin1,spin2,l,im,uu,ud,du,dd) &
156 6 : !$OMP private(m,mp) collapse(2)
157 : DO m = -l,l
158 : DO mp = -l,l
159 : im(:,m,mp,jspin) = uu(:,m,mp,jspin) * scalarGF%uun(spin1,spin2) &
160 : + ud(:,m,mp,jspin) * scalarGF%udn(spin1,spin2) &
161 : + du(:,m,mp,jspin) * scalarGF%dun(spin1,spin2) &
162 : + dd(:,m,mp,jspin) * scalarGF%ddn(spin1,spin2)
163 : ENDDO
164 : ENDDO
165 : !$OMP end parallel do
166 : ENDDO
167 :
168 2 : CALL kk_cutoff(im,noco,l_mperp,l,input%jspins,eMesh,cutoff,scalingFactor)
169 :
170 :
171 2 : END SUBROUTINE kk_cutoffRadial
172 :
173 2 : SUBROUTINE kk_cutoffRadialLO(uu,ud,du,dd,uulo,dulo,ulou,ulod,uloulop,atoms,noco,scalarGF,l_mperp,&
174 2 : l,atomtype,input,eMesh,cutoff,scalingFactor)
175 :
176 : COMPLEX, INTENT(IN) :: uu(:,-lmaxU_const:,-lmaxU_const:,:)
177 : COMPLEX, INTENT(IN) :: ud(:,-lmaxU_const:,-lmaxU_const:,:)
178 : COMPLEX, INTENT(IN) :: du(:,-lmaxU_const:,-lmaxU_const:,:)
179 : COMPLEX, INTENT(IN) :: dd(:,-lmaxU_const:,-lmaxU_const:,:)
180 : COMPLEX, INTENT(IN) :: uulo(:,-lmaxU_const:,-lmaxU_const:,:,:)
181 : COMPLEX, INTENT(IN) :: dulo(:,-lmaxU_const:,-lmaxU_const:,:,:)
182 : COMPLEX, INTENT(IN) :: ulou(:,-lmaxU_const:,-lmaxU_const:,:,:)
183 : COMPLEX, INTENT(IN) :: ulod(:,-lmaxU_const:,-lmaxU_const:,:,:)
184 : COMPLEX, INTENT(IN) :: uloulop(:,-lmaxU_const:,-lmaxU_const:,:,:,:)
185 : TYPE(t_noco), INTENT(IN) :: noco
186 : type(t_atoms), intent(in) :: atoms
187 : TYPE(t_scalarGF), INTENT(IN) :: scalarGF
188 : LOGICAL, INTENT(IN) :: l_mperp
189 : INTEGER, INTENT(IN) :: l, atomtype
190 : TYPE(t_input), INTENT(IN) :: input
191 : REAL, INTENT(IN) :: eMesh(:)
192 : INTEGER, INTENT(INOUT) :: cutoff(:,:)
193 : REAL, INTENT(INOUT) :: scalingFactor(:)
194 :
195 2 : COMPLEX, ALLOCATABLE :: im(:,:,:,:)
196 :
197 : INTEGER :: jspin,m,mp,spin1,spin2,ilo,ilop,iLO_ind,iLOp_ind
198 :
199 : !calculate the spherical average from the original greens function
200 1058638 : ALLOCATE(im(SIZE(uu,1),-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const,SIZE(uu,4)),source=cmplx_0)
201 6 : DO jspin = 1, SIZE(im,4)
202 4 : IF(jspin < 3) THEN
203 4 : spin1 = jspin
204 4 : spin2 = jspin
205 : ELSE
206 : spin1 = 2
207 : spin2 = 1
208 : ENDIF
209 : !$OMP parallel do default(none) &
210 : !$OMP shared(scalarGF,jspin,spin1,spin2,l,im,uu,ud,du,dd,atoms) &
211 : !$OMP shared(uulo,ulou,ulod,dulo,uloulop,atomtype) &
212 6 : !$OMP private(m,mp,ilo,ilop,iLO_ind,iLOp_ind) collapse(2)
213 : DO m = -l,l
214 : DO mp = -l,l
215 : im(:,m,mp,jspin) = uu(:,m,mp,jspin) * scalarGF%uun(spin1,spin2) &
216 : + ud(:,m,mp,jspin) * scalarGF%udn(spin1,spin2) &
217 : + du(:,m,mp,jspin) * scalarGF%dun(spin1,spin2) &
218 : + dd(:,m,mp,jspin) * scalarGF%ddn(spin1,spin2)
219 : iLO_ind = 0
220 : DO ilo = 1, atoms%nlo(atomType)
221 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
222 : iLO_ind = iLO_ind + 1
223 : im(:,m,mp,jspin) = im(:,m,mp,jspin) &
224 : + uulo(:,m,mp,iLO_ind,jspin) * scalarGF%uulon(ilo,spin1,spin2) &
225 : + dulo(:,m,mp,iLO_ind,jspin) * scalarGF%dulon(ilo,spin1,spin2) &
226 : + ulou(:,m,mp,iLO_ind,jspin) * scalarGF%uloun(ilo,spin1,spin2) &
227 : + ulod(:,m,mp,iLO_ind,jspin) * scalarGF%ulodn(ilo,spin1,spin2)
228 : iLOp_ind = 0
229 : DO ilop = 1, atoms%nlo(atomType)
230 : IF(atoms%llo(ilop,atomType).NE.l) CYCLE
231 : iLOp_ind = iLOp_ind + 1
232 : im(:,m,mp,jspin) = im(:,m,mp,jspin) &
233 : + uloulop(:,m,mp,iLO_ind,iLOp_ind,jspin) * scalarGF%uloulopn(ilo,ilop,spin1,spin2)
234 : ENDDO
235 : enddo
236 : enddo
237 : ENDDO
238 : !$OMP end parallel do
239 : ENDDO
240 :
241 2 : CALL kk_cutoff(im,noco,l_mperp,l,input%jspins,eMesh,cutoff,scalingFactor)
242 :
243 2 : END SUBROUTINE kk_cutoffRadialLO
244 :
245 : END MODULE m_kk_cutoff
|
