Line data Source code
1 : MODULE m_greensfCalcRealPart
2 :
3 : !------------------------------------------------------------------------------
4 : !
5 : ! MODULE: m_greensfCalcRealPart
6 : !
7 : !> @author
8 : !> Henning Janßen
9 : !
10 : ! DESCRIPTION:
11 : !> This module contains the functions to calculate the imaginary part of the
12 : !> onsite GF with and without radial dependence
13 : !> Further we can transform this imaginary part to obtain the Green's Function
14 : !> using the Kramer Kronig Transformation
15 : !
16 : !------------------------------------------------------------------------------
17 :
18 : USE m_juDFT
19 : USE m_types
20 : USE m_constants
21 : USE m_kkintgr
22 : USE m_kk_cutoff
23 :
24 : IMPLICIT NONE
25 :
26 : private
27 : public greensfCalcRealPart
28 :
29 : CONTAINS
30 :
31 42 : SUBROUTINE greensfCalcRealPart(atoms,gfinp,sym,input,enpara,noco,kpts,fmpi,ef,greensfImagPart,g)
32 :
33 : TYPE(t_atoms), INTENT(IN) :: atoms
34 : TYPE(t_gfinp), INTENT(IN) :: gfinp
35 : TYPE(t_sym), INTENT(IN) :: sym
36 : TYPE(t_noco), INTENT(IN) :: noco
37 : TYPE(t_kpts), INTENT(IN) :: kpts
38 : TYPE(t_input), INTENT(IN) :: input
39 : TYPE(t_mpi), INTENT(IN) :: fmpi
40 : type(t_enpara), intent(in) :: enpara
41 : REAL, INTENT(IN) :: ef
42 : TYPE(t_greensfImagPart), INTENT(INOUT) :: greensfImagPart
43 : TYPE(t_greensf), INTENT(INOUT) :: g(:)
44 :
45 : INTEGER :: i_gf,i_elem,l,m,mp,indUnique,nLO,iLO,iLOp,i_elemLO
46 : INTEGER :: jspin,nspins,ipm,lp,refCutoff
47 : INTEGER :: contourShape, iContour
48 : INTEGER :: i_gf_start,i_gf_end,spin_start,spin_end
49 : INTEGER :: ikpt, ikpt_i, refElem
50 : LOGICAL :: l_fixedCutoffset,l_sphavg,l_kresolved_int,l_kresolved
51 : REAL :: del,eb,fixedCutoff,bk(3)
52 42 : REAL, ALLOCATABLE :: eMesh(:)
53 42 : COMPLEX, ALLOCATABLE :: gmat(:,:,:),imag(:,:,:)
54 :
55 : !Get the information on the real axis energy mesh
56 42 : CALL gfinp%eMesh(ef,del,eb,eMesh=eMesh)
57 :
58 42 : nspins = MERGE(3,input%jspins,gfinp%l_mperp)
59 :
60 42 : IF(fmpi%irank.EQ.0) THEN
61 21 : CALL timestart("Green's Function: Integration Cutoff")
62 530 : DO i_gf = 1, gfinp%n
63 :
64 : !Get the information of ith current element
65 509 : l = g(i_gf)%elem%l
66 509 : lp = g(i_gf)%elem%lp
67 509 : l_sphavg = g(i_gf)%elem%l_sphavg
68 509 : l_fixedCutoffset = g(i_gf)%elem%l_fixedCutoffset
69 509 : fixedCutoff = g(i_gf)%elem%fixedCutoff
70 509 : refCutoff = g(i_gf)%elem%refCutoff
71 509 : l_kresolved_int = g(i_gf)%elem%l_kresolved_int
72 :
73 509 : IF(refCutoff /= -1) CYCLE
74 :
75 31 : IF(l_fixedCutoffset) THEN
76 0 : greensfImagPart%kkintgr_cutoff(i_gf,:,1) = 1
77 0 : greensfImagPart%kkintgr_cutoff(i_gf,:,2) = INT((fixedCutoff+ef-eb)/del)+1
78 : CYCLE
79 : ENDIF
80 :
81 52 : IF(.NOT.gfinp%isUnique(i_gf,distinct_kresolved_int=.TRUE.)) THEN
82 2 : indUnique = gfinp%getUniqueElement(i_gf,distinct_kresolved_int=.TRUE.)
83 : !This cutoff was already calculated
84 14 : greensfImagPart%kkintgr_cutoff(i_gf,:,:) = greensfImagPart%kkintgr_cutoff(indUnique,:,:)
85 : ELSE
86 29 : i_elem = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg,l_kresolved_int=l_kresolved_int)
87 : IF(.not. g(i_gf)%elem%isOffDiag() .and. &
88 29 : .not. has_sclos(g(i_gf)%elem, atoms,enpara) .and. &
89 478 : .not. l_kresolved_int) THEN
90 : !
91 : !Check the integral over the PDOS to define a cutoff for the Kramer-Kronigs-Integration
92 : ! with SCLOs I just use a fixed cutoff or reference otherwise I would need to check whether
93 : ! the SCLO lies in the energy boundary and raise the expected number of states accordingly
94 26 : IF(l_sphavg) THEN
95 : CALL kk_cutoff(greensfImagPart%sphavg(:,:,:,i_elem,:),noco,gfinp%l_mperp,l,input%jspins,&
96 22 : eMesh,greensfImagPart%kkintgr_cutoff(i_gf,:,:),greensfImagPart%scalingFactorSphavg(i_elem,:))
97 4 : ELSE IF(g(i_gf)%elem%countLOs(atoms)==0) then
98 : !Onsite element with radial dependence
99 : CALL kk_cutoffRadial(greensfImagPart%uu(:,:,:,i_elem,:),greensfImagPart%ud(:,:,:,i_elem,:),&
100 : greensfImagPart%du(:,:,:,i_elem,:),greensfImagPart%dd(:,:,:,i_elem,:),&
101 : noco,g(i_gf)%scalarProducts,gfinp%l_mperp,l,input,eMesh,&
102 2 : greensfImagPart%kkintgr_cutoff(i_gf,:,:),greensfImagPart%scalingFactorRadial(i_elem,:))
103 : ELSE
104 : CALL kk_cutoffRadialLO(greensfImagPart%uu(:,:,:,i_elem,:),greensfImagPart%ud(:,:,:,i_elem,:),&
105 : greensfImagPart%du(:,:,:,i_elem,:),greensfImagPart%dd(:,:,:,i_elem,:),&
106 : greensfImagPart%uulo(:,:,:,:,i_elem,:),greensfImagPart%dulo(:,:,:,:,i_elem,:),&
107 : greensfImagPart%ulou(:,:,:,:,i_elem,:),greensfImagPart%ulod(:,:,:,:,i_elem,:),&
108 : greensfImagPart%uloulop(:,:,:,:,:,i_elem,:),&
109 : atoms,noco,g(i_gf)%scalarProducts,gfinp%l_mperp,l,g(i_gf)%elem%atomType,input,eMesh,&
110 2 : greensfImagPart%kkintgr_cutoff(i_gf,:,:),greensfImagPart%scalingFactorRadial(i_elem,:))
111 : ENDIF
112 : ELSE
113 : !For all other elements we just use ef+elup as a hard cutoff
114 9 : greensfImagPart%kkintgr_cutoff(i_gf,:,1) = 1
115 9 : greensfImagPart%kkintgr_cutoff(i_gf,:,2) = gfinp%ne
116 : ENDIF
117 : ENDIF
118 :
119 : ENDDO
120 :
121 : !Getting reference Cutoffs and perform scaling
122 530 : DO i_gf = 1, gfinp%n
123 509 : refCutoff = g(i_gf)%elem%refCutoff
124 509 : l_kresolved_int = g(i_gf)%elem%l_kresolved_int
125 509 : l_sphavg = g(i_gf)%elem%l_sphavg
126 509 : i_elem = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg,l_kresolved_int=l_kresolved_int)
127 509 : i_elemLO = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg,l_kresolved_int=l_kresolved_int,lo=.TRUE.)
128 509 : nLO = g(i_gf)%elem%countLOs(atoms)
129 :
130 509 : IF(refCutoff/=-1) THEN
131 : !Overwrite cutoff with reference from other elements
132 3346 : greensfImagPart%kkintgr_cutoff(i_gf,:,:) = greensfImagPart%kkintgr_cutoff(refCutoff,:,:)
133 :
134 : ! refElem = gfinp%getUniqueElement(refCutoff,distinct_kresolved_int=.TRUE.)
135 : ! if (l_sphavg .and. .not. l_kresolved_int) then
136 : ! greensfImagPart%scalingFactorSphavg(i_elem,:) = greensfImagPart%scalingFactorSphavg(refElem,:)
137 : ! else if(.not.l_sphavg) then
138 : ! greensfImagPart%scalingFactorRadial(i_elem,:) = greensfImagPart%scalingFactorRadial(refElem,:)
139 : ! else
140 : ! greensfImagPart%scalingFactorSphavgKres(i_elem,:) = greensfImagPart%scalingFactorSphavgKres(refElem,:)
141 : ! endif
142 : ENDIF
143 509 : if(.NOT.gfinp%isUnique(i_gf,distinct_kresolved_int=.TRUE.)) cycle
144 151 : CALL greensfImagPart%scale(i_elem,i_elemLO,l_sphavg,nLO,k_resolved=l_kresolved_int)
145 : ENDDO
146 21 : CALL timestop("Green's Function: Integration Cutoff")
147 : ENDIF
148 :
149 : !Broadcast cutoffs and modified imaginary parts
150 42 : CALL greensfImagPart%mpi_bc(fmpi%mpi_comm)
151 :
152 : !Distribute the Calculations
153 42 : CALL gfinp%distribute_elements(fmpi%irank, fmpi%isize, nspins, i_gf_start, i_gf_end, spin_start, spin_end)
154 :
155 : !Initialize kkintgr_module variables
156 559 : DO i_gf = i_gf_start, i_gf_end
157 517 : IF(i_gf.LT.1 .OR. i_gf.GT.gfinp%n) CYCLE !Make sure to not produce segfaults with mpi
158 517 : contourShape = gfinp%contour(g(i_gf)%elem%iContour)%shape
159 :
160 559 : CALL kkintgr_init(eMesh,g(i_gf)%contour%e,g(i_gf)%elem%iContour,gfinp%numberContours, contourShape, gfinp%additional_smearing)
161 : ENDDO
162 :
163 559 : DO i_gf = i_gf_start, i_gf_end
164 :
165 517 : IF(i_gf.LT.1 .OR. i_gf.GT.gfinp%n) CYCLE !Make sure to not produce segfaults with mpi
166 :
167 : !Get the information of ith current element
168 517 : l = g(i_gf)%elem%l
169 517 : lp = g(i_gf)%elem%lp
170 517 : l_sphavg = g(i_gf)%elem%l_sphavg
171 517 : iContour = g(i_gf)%elem%iContour
172 517 : nLO = g(i_gf)%elem%countLOs(atoms)
173 517 : IF(g(i_gf)%elem%representative_elem > 0) CYCLE
174 140 : IF(g(i_gf)%elem%l_kresolved_int) CYCLE
175 :
176 140 : i_elem = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg,l_kresolved_int=.FALSE.)
177 140 : i_elemLO = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg,lo=.TRUE.,l_kresolved_int=.FALSE.)
178 :
179 140 : CALL timestart("Green's Function: Kramer-Kronigs-Integration")
180 408 : DO jspin = spin_start, spin_end
181 944 : DO ipm = 1, 2 !upper or lower half of the complex plane (G(E \pm i delta))
182 :
183 804 : IF(l_sphavg) THEN
184 183897296 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,jspin,l_sphavg)
185 512 : CALL kkintgr(imag,ipm==2,g(i_gf)%gmmpMat(:,:,:,jspin,ipm),iContour)
186 : ELSE
187 : ! In the case of radial dependence we perform the kramers-kronig-integration seperately for uu,dd,etc.
188 : ! We can do this because the radial functions are independent of E
189 6351792 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,jspin,l_sphavg,imat=1)
190 24 : CALL kkintgr(imag,ipm==2,g(i_gf)%uu(:,:,:,jspin,ipm),iContour)
191 6351792 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,jspin,l_sphavg,imat=2)
192 24 : CALL kkintgr(imag,ipm==2,g(i_gf)%dd(:,:,:,jspin,ipm),iContour)
193 6351792 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,jspin,l_sphavg,imat=3)
194 24 : CALL kkintgr(imag,ipm==2,g(i_gf)%ud(:,:,:,jspin,ipm),iContour)
195 6351792 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,jspin,l_sphavg,imat=4)
196 24 : CALL kkintgr(imag,ipm==2,g(i_gf)%du(:,:,:,jspin,ipm),iContour)
197 :
198 : !KKT for LOs
199 24 : IF(nLO>0) THEN
200 36 : DO iLO = 1, nLO
201 5293160 : imag = greensfImagPart%applyCutoff(i_elemLO,i_gf,jspin,l_sphavg,imat=1,iLO=iLO)
202 20 : CALL kkintgr(imag,ipm==2,g(i_gf)%uulo(:,:,:,iLO,jspin,ipm),iContour)
203 5293160 : imag = greensfImagPart%applyCutoff(i_elemLO,i_gf,jspin,l_sphavg,imat=2,iLO=iLO)
204 20 : CALL kkintgr(imag,ipm==2,g(i_gf)%ulou(:,:,:,iLO,jspin,ipm),iContour)
205 5293160 : imag = greensfImagPart%applyCutoff(i_elemLO,i_gf,jspin,l_sphavg,imat=3,iLO=iLO)
206 20 : CALL kkintgr(imag,ipm==2,g(i_gf)%dulo(:,:,:,iLO,jspin,ipm),iContour)
207 5293160 : imag = greensfImagPart%applyCutoff(i_elemLO,i_gf,jspin,l_sphavg,imat=4,iLO=iLO)
208 20 : CALL kkintgr(imag,ipm==2,g(i_gf)%ulod(:,:,:,iLO,jspin,ipm),iContour)
209 :
210 64 : DO iLOp = 1, nLO
211 7410424 : imag = greensfImagPart%applyCutoff(i_elemLO,i_gf,jspin,l_sphavg,iLO=iLO,iLOp=iLop)
212 48 : CALL kkintgr(imag,ipm==2,g(i_gf)%uloulop(:,:,:,iLO,iLOp,jspin,ipm),iContour)
213 : ENDDO
214 : ENDDO
215 : ENDIF
216 : ENDIF
217 : ENDDO
218 : ENDDO
219 559 : CALL timestop("Green's Function: Kramer-Kronigs-Integration")
220 : ENDDO
221 42 : CALL kkintgr_free()
222 :
223 1060 : IF(ANY(gfinp%elem(:)%l_kresolved_int)) THEN
224 0 : CALL gfinp%distribute_elements(fmpi%n_rank, fmpi%n_size, nspins, i_gf_start, i_gf_end, spin_start, spin_end, k_resolved=.TRUE.)
225 : !Initialize kkintgr_module variables
226 0 : DO i_gf = i_gf_start, i_gf_end
227 0 : IF(i_gf.LT.1 .OR. i_gf.GT.gfinp%n) CYCLE !Make sure to not produce segfaults with mpi
228 0 : contourShape = gfinp%contour(g(i_gf)%elem%iContour)%shape
229 0 : CALL kkintgr_init(eMesh,g(i_gf)%contour%e,g(i_gf)%elem%iContour,gfinp%numberContours, contourShape, gfinp%additional_smearing)
230 : ENDDO
231 0 : CALL timestart("Green's Function: K-Resolved Kramer-Kronigs-Integration")
232 0 : DO ikpt_i = 1, SIZE(fmpi%k_list)
233 0 : ikpt = fmpi%k_list(ikpt_i)
234 0 : bk = kpts%bk(:,ikpt)
235 0 : DO i_gf = i_gf_start, i_gf_end
236 :
237 0 : IF(i_gf.LT.1 .OR. i_gf.GT.gfinp%n) CYCLE !Make sure to not produce segfaults with mpi
238 :
239 : !Get the information of ith current element
240 0 : l = g(i_gf)%elem%l
241 0 : lp = g(i_gf)%elem%lp
242 0 : l_sphavg = g(i_gf)%elem%l_sphavg
243 0 : l_kresolved = g(i_gf)%elem%l_kresolved
244 0 : iContour = g(i_gf)%elem%iContour
245 0 : nLO = g(i_gf)%elem%countLOs(atoms)
246 0 : IF(g(i_gf)%elem%representative_elem > 0) CYCLE
247 0 : IF(.NOT.g(i_gf)%elem%l_kresolved_int) CYCLE
248 :
249 0 : i_elem = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg,l_kresolved_int=.TRUE.)
250 0 : i_elemLO = gfinp%uniqueElements(atoms,max_index=i_gf,l_sphavg=l_sphavg,lo=.TRUE.,l_kresolved_int=.TRUE.)
251 :
252 0 : IF(.NOT.l_kresolved) ALLOCATE(gmat(SIZE(g(i_gf)%contour%e),-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const), source=cmplx_0)
253 :
254 0 : CALL timestart("Green's Function: Kramer-Kronigs-Integration")
255 0 : DO jspin = spin_start, spin_end
256 0 : DO ipm = 1, 2 !upper or lower half of the complex plane (G(E \pm i delta))
257 :
258 0 : IF(l_sphavg) THEN
259 0 : imag = greensfImagPart%applyCutoff(i_elem,i_gf,jspin,l_sphavg,ikpt=ikpt_i)
260 0 : IF(l_kresolved) THEN
261 0 : CALL kkintgr(imag,ipm==2,g(i_gf)%gmmpMat_k(:,:,:,jspin,ipm,ikpt),iContour)
262 : ELSE
263 0 : CALL kkintgr(imag,ipm==2,gmat,iContour)
264 0 : g(i_gf)%gmmpMat(:,:,:,jspin,ipm) = g(i_gf)%gmmpMat(:,:,:,jspin,ipm) + gmat
265 : ENDIF
266 : ELSE
267 0 : CALL juDFT_error("No Green's function with k-resolution and radial dependence implemented")
268 : ENDIF
269 :
270 : ENDDO
271 : ENDDO
272 0 : CALL timestop("Green's Function: Kramer-Kronigs-Integration")
273 0 : IF(.NOT.l_kresolved) DEALLOCATE(gmat)
274 : ENDDO
275 : ENDDO
276 0 : CALL timestop("Green's Function: K-Resolved Kramer-Kronigs-Integration")
277 0 : CALL kkintgr_free()
278 : ENDIF
279 :
280 :
281 : #ifdef CPP_MPI
282 42 : CALL timestart("Green's Function: Collect")
283 : !Collect all the greensFuntions
284 1060 : DO i_gf = 1, gfinp%n
285 1060 : CALL g(i_gf)%collect(fmpi%mpi_comm)
286 : ENDDO
287 42 : CALL timestop("Green's Function: Collect")
288 : #endif
289 :
290 42 : IF(fmpi%irank.EQ.0) THEN
291 : !perform rotations for intersite elements
292 530 : DO i_gf = 1, gfinp%n
293 509 : IF(g(i_gf)%elem%representative_elem <= 0) CYCLE
294 377 : CALL g(i_gf)%set_gfdata(g(g(i_gf)%elem%representative_elem))
295 530 : CALL g(i_gf)%rotate(sym,atoms)
296 : ENDDO
297 : ENDIF
298 :
299 1060 : DO i_gf = 1, gfinp%n
300 1060 : CALL g(i_gf)%mpi_bc(fmpi%mpi_comm)
301 : ENDDO
302 :
303 42 : END SUBROUTINE greensfCalcRealPart
304 :
305 29 : logical function has_sclos(elem,atoms, enpara)
306 :
307 : type(t_gfelementtype), intent(in) :: elem
308 : type(t_atoms), intent(in) :: atoms
309 : type(t_enpara), intent(in) :: enpara
310 :
311 : integer :: ilo
312 :
313 29 : has_sclos = .false.
314 93 : DO ilo = 1, atoms%nlo(elem%atomType)
315 64 : if(atoms%llo(ilo,elem%atomType).NE.elem%l) cycle
316 : !Check for HELO (negative energy parameter) and HDLO (eDeriv > 0)
317 11 : if(all(enpara%qn_ello(ilo,elem%atomType,:)<0).or.atoms%ulo_der(ilo, elem%atomType)>=1) cycle
318 154 : has_sclos = .true.
319 : ENDDO
320 :
321 29 : IF(elem%l.NE.elem%lp.OR.elem%atomType.NE.elem%atomTypep) THEN
322 0 : DO ilo = 1, atoms%nlo(elem%atomTypep)
323 0 : if(atoms%llo(ilo,elem%atomTypep).NE.elem%lp) cycle
324 : !Check for HELO (negative energy parameter) and HDLO (eDeriv > 0)
325 0 : if(all(enpara%qn_ello(ilo,elem%atomTypep,:)<0).or.atoms%ulo_der(ilo, elem%atomTypep)>=1) cycle
326 29 : has_sclos = .true.
327 : ENDDO
328 : ENDIF
329 :
330 29 : end function
331 :
332 : END MODULE m_greensfCalcRealPart
|
