Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : MODULE m_types_greensf
8 :
9 : !------------------------------------------------------------------------------
10 : !
11 : ! MODULE: m_types_greensf
12 : !
13 : !> @author
14 : !> Henning Janßen
15 : !
16 : ! DESCRIPTION:
17 : !> Contains a type for onsite and intersite green's functions in the mt-sphere
18 : !> It stores the energy contour in the complex plane and the corresponding
19 : !> matrix elements of the green's function
20 : !------------------------------------------------------------------------------
21 :
22 : USE m_juDFT
23 : USE m_constants
24 : USE m_types_setup
25 : USE m_types_greensfContourData
26 : USE m_types_scalarGF
27 : USE m_types_nococonv
28 :
29 : IMPLICIT NONE
30 :
31 : PRIVATE
32 :
33 : TYPE t_greensf
34 :
35 : LOGICAL :: l_calc = .FALSE.
36 : LOGICAL :: l_sphavg
37 : LOGICAL :: l_kresolved
38 :
39 : !Energy contour parameters
40 : TYPE(t_greensfContourData) :: contour
41 :
42 : !Pointer to the element type in gfinp
43 : TYPE(t_gfelementtype), POINTER :: elem => NULL()
44 : TYPE(t_scalarGF) :: scalarProducts
45 :
46 : !Arrays for Green's function
47 : COMPLEX, ALLOCATABLE :: gmmpMat(:,:,:,:,:)
48 :
49 : !for radial dependence
50 : COMPLEX, ALLOCATABLE :: uu(:,:,:,:,:)
51 : COMPLEX, ALLOCATABLE :: dd(:,:,:,:,:)
52 : COMPLEX, ALLOCATABLE :: du(:,:,:,:,:)
53 : COMPLEX, ALLOCATABLE :: ud(:,:,:,:,:)
54 :
55 : COMPLEX, ALLOCATABLE :: uulo(:,:,:,:,:,:)
56 : COMPLEX, ALLOCATABLE :: ulou(:,:,:,:,:,:)
57 : COMPLEX, ALLOCATABLE :: dulo(:,:,:,:,:,:)
58 : COMPLEX, ALLOCATABLE :: ulod(:,:,:,:,:,:)
59 :
60 : COMPLEX, ALLOCATABLE :: uloulop(:,:,:,:,:,:,:)
61 :
62 : COMPLEX, ALLOCATABLE :: gmmpMat_k(:,:,:,:,:,:)
63 :
64 : CONTAINS
65 : PROCEDURE, PASS :: init => init_greensf
66 : PROCEDURE :: mpi_bc => mpi_bc_greensf
67 : PROCEDURE :: collect => collect_greensf
68 : PROCEDURE :: get => get_gf
69 : PROCEDURE :: occmtx_greensf_spin
70 : PROCEDURE :: occmtx_greensf_complete
71 : GENERIC :: occupationMatrix => occmtx_greensf_spin, occmtx_greensf_complete
72 : PROCEDURE :: getFullMatrix => getFullMatrix_gf
73 : PROCEDURE :: getRadial => getRadial_gf
74 : PROCEDURE :: getRadialRadial => getRadialRadial_gf!(Full Radial dependence for intersite)
75 : PROCEDURE :: integrateOverMT => integrateOverMT_greensf
76 : PROCEDURE :: set => set_gf
77 : PROCEDURE :: set_gfdata => set_gfdata
78 : PROCEDURE :: rotate => rotate_gf
79 : PROCEDURE :: rotate_euler_angles => rotate_euler_angles_gf
80 : PROCEDURE :: reset => reset_gf
81 : PROCEDURE :: resetSingleElem => resetSingleElem_gf
82 : PROCEDURE :: checkEmpty => checkEmpty_greensf
83 : END TYPE t_greensf
84 :
85 : PUBLIC t_greensf
86 :
87 : CONTAINS
88 :
89 1012 : SUBROUTINE init_greensf(this,gfelem,gfinp,atoms,input,contour_in,l_sphavg_in,l_kresolved_in,nkpt)
90 :
91 : CLASS(t_greensf), INTENT(INOUT) :: this
92 : TYPE(t_gfelementtype), TARGET, INTENT(IN) :: gfelem
93 : TYPE(t_gfinp), INTENT(IN) :: gfinp
94 : TYPE(t_atoms), INTENT(IN) :: atoms
95 : TYPE(t_input), INTENT(IN) :: input
96 : !Pass a already calculated energy contour to the type
97 : TYPE(t_greensfContourData), OPTIONAL, INTENT(IN) :: contour_in
98 : LOGICAL, OPTIONAL, INTENT(IN) :: l_sphavg_in !To overwrite the allocation for integrateOverMT_greensf
99 : LOGICAL, OPTIONAL, INTENT(IN) :: l_kresolved_in !To overwrite the allocation for integrateOverBZ_greensf
100 : INTEGER, OPTIONAL, INTENT(IN) :: nkpt
101 :
102 : INTEGER spin_dim,lmax,nLO
103 : LOGICAL l_sphavg
104 :
105 1012 : this%elem => gfelem
106 1012 : this%l_calc = .FALSE.
107 :
108 1012 : nLO = this%elem%countLOs(atoms)
109 :
110 : !Initialize the contour
111 1012 : CALL this%contour%init(gfinp%contour(this%elem%iContour),contour_in=contour_in)
112 :
113 1012 : spin_dim = MERGE(3,input%jspins,gfinp%l_mperp)
114 1012 : lmax = lmaxU_const
115 :
116 1012 : this%l_sphavg = this%elem%l_sphavg
117 1012 : IF(PRESENT(l_sphavg_in)) this%l_sphavg = l_sphavg_in
118 :
119 1012 : this%l_kresolved = this%elem%l_kresolved
120 1012 : IF(PRESENT(l_kresolved_in)) this%l_kresolved = l_kresolved_in
121 :
122 1012 : IF(this%l_kresolved.AND. .NOT.PRESENT(nkpt)) CALL juDFT_error("Missing argument nkpt for k-resolved green's function", calledby='init_greensf')
123 :
124 1012 : IF(this%l_sphavg) THEN
125 1000 : IF(this%l_kresolved) THEN
126 0 : ALLOCATE(this%gmmpMat_k(this%contour%nz,-lmax:lmax,-lmax:lmax,spin_dim,2,nkpt),source=cmplx_0)
127 : ELSE
128 27509528 : ALLOCATE(this%gmmpMat(this%contour%nz,-lmax:lmax,-lmax:lmax,spin_dim,2),source=cmplx_0)
129 : ENDIF
130 : ELSE
131 303888 : ALLOCATE(this%uu(this%contour%nz,-lmax:lmax,-lmax:lmax,spin_dim,2),source=cmplx_0)
132 303852 : ALLOCATE(this%dd(this%contour%nz,-lmax:lmax,-lmax:lmax,spin_dim,2),source=cmplx_0)
133 303852 : ALLOCATE(this%du(this%contour%nz,-lmax:lmax,-lmax:lmax,spin_dim,2),source=cmplx_0)
134 303852 : ALLOCATE(this%ud(this%contour%nz,-lmax:lmax,-lmax:lmax,spin_dim,2),source=cmplx_0)
135 :
136 12 : IF(nLO>0) THEN
137 253256 : ALLOCATE(this%uulo(this%contour%nz,-lmax:lmax,-lmax:lmax,nLO,spin_dim,2),source=cmplx_0)
138 253240 : ALLOCATE(this%ulou(this%contour%nz,-lmax:lmax,-lmax:lmax,nLO,spin_dim,2),source=cmplx_0)
139 253240 : ALLOCATE(this%dulo(this%contour%nz,-lmax:lmax,-lmax:lmax,nLO,spin_dim,2),source=cmplx_0)
140 253240 : ALLOCATE(this%ulod(this%contour%nz,-lmax:lmax,-lmax:lmax,nLO,spin_dim,2),source=cmplx_0)
141 :
142 354568 : ALLOCATE(this%uloulop(this%contour%nz,-lmax:lmax,-lmax:lmax,nLO,nLO,spin_dim,2),source=cmplx_0)
143 : ENDIF
144 : ENDIF
145 :
146 1012 : END SUBROUTINE init_greensf
147 :
148 1018 : SUBROUTINE mpi_bc_greensf(this,mpi_comm,irank)
149 : USE m_mpi_bc_tool
150 : CLASS(t_greensf), INTENT(INOUT)::this
151 : INTEGER, INTENT(IN):: mpi_comm
152 : INTEGER, INTENT(IN), OPTIONAL::irank
153 : INTEGER ::rank
154 1018 : IF (PRESENT(irank)) THEN
155 0 : rank = irank
156 : ELSE
157 1018 : rank = 0
158 : END IF
159 :
160 1018 : CALL mpi_bc(this%l_calc,rank,mpi_comm)
161 1018 : CALL mpi_bc(this%l_sphavg,rank,mpi_comm)
162 1018 : CALL mpi_bc(this%l_kresolved,rank,mpi_comm)
163 :
164 1018 : CALL this%contour%mpi_bc(mpi_comm,rank)
165 :
166 1018 : IF(ALLOCATED(this%gmmpMat)) CALL mpi_bc(this%gmmpMat,rank,mpi_comm)
167 1018 : IF(ALLOCATED(this%gmmpMat_k)) CALL mpi_bc(this%gmmpMat_k,rank,mpi_comm)
168 1018 : IF(ALLOCATED(this%uu)) CALL mpi_bc(this%uu,rank,mpi_comm)
169 1018 : IF(ALLOCATED(this%ud)) CALL mpi_bc(this%ud,rank,mpi_comm)
170 1018 : IF(ALLOCATED(this%du)) CALL mpi_bc(this%du,rank,mpi_comm)
171 1018 : IF(ALLOCATED(this%dd)) CALL mpi_bc(this%dd,rank,mpi_comm)
172 1018 : IF(ALLOCATED(this%uulo)) CALL mpi_bc(this%uulo,rank,mpi_comm)
173 1018 : IF(ALLOCATED(this%ulou)) CALL mpi_bc(this%ulou,rank,mpi_comm)
174 1018 : IF(ALLOCATED(this%dulo)) CALL mpi_bc(this%dulo,rank,mpi_comm)
175 1018 : IF(ALLOCATED(this%ulod)) CALL mpi_bc(this%ulod,rank,mpi_comm)
176 1018 : IF(ALLOCATED(this%uloulop)) CALL mpi_bc(this%uloulop,rank,mpi_comm)
177 :
178 1018 : END SUBROUTINE mpi_bc_greensf
179 :
180 1018 : SUBROUTINE collect_greensf(this,mpi_communicator)
181 :
182 : #ifdef CPP_MPI
183 : USE mpi
184 : #endif
185 :
186 : CLASS(t_greensf), INTENT(INOUT) :: this
187 : INTEGER, INTENT(IN) :: mpi_communicator
188 : #ifdef CPP_MPI
189 : INTEGER:: ierr,n
190 1018 : COMPLEX,ALLOCATABLE::ctmp(:)
191 :
192 1018 : IF(ALLOCATED(this%gmmpMat)) THEN
193 6036 : n = SIZE(this%gmmpMat)
194 3018 : ALLOCATE(ctmp(n))
195 1006 : CALL MPI_ALLREDUCE(this%gmmpMat,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
196 1006 : CALL zcopy(n,ctmp,1,this%gmmpMat,1)
197 1006 : DEALLOCATE(ctmp)
198 12 : ELSE IF(ALLOCATED(this%gmmpMat_k)) THEN
199 0 : n = SIZE(this%gmmpMat_k)
200 0 : ALLOCATE(ctmp(n))
201 0 : CALL MPI_ALLREDUCE(this%gmmpMat_k,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
202 0 : CALL zcopy(n,ctmp,1,this%gmmpMat_k,1)
203 0 : DEALLOCATE(ctmp)
204 : ELSE
205 72 : n = SIZE(this%uu)
206 36 : ALLOCATE(ctmp(n))
207 12 : CALL MPI_ALLREDUCE(this%uu,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
208 12 : CALL zcopy(n,ctmp,1,this%uu,1)
209 12 : CALL MPI_ALLREDUCE(this%ud,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
210 12 : CALL zcopy(n,ctmp,1,this%ud,1)
211 12 : CALL MPI_ALLREDUCE(this%du,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
212 12 : CALL zcopy(n,ctmp,1,this%du,1)
213 12 : CALL MPI_ALLREDUCE(this%dd,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
214 12 : CALL zcopy(n,ctmp,1,this%dd,1)
215 12 : DEALLOCATE(ctmp)
216 :
217 12 : IF(ALLOCATED(this%uulo)) THEN
218 56 : n = SIZE(this%uulo)
219 24 : ALLOCATE(ctmp(n))
220 8 : CALL MPI_ALLREDUCE(this%uulo,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
221 8 : CALL zcopy(n,ctmp,1,this%uulo,1)
222 8 : CALL MPI_ALLREDUCE(this%ulou,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
223 8 : CALL zcopy(n,ctmp,1,this%ulou,1)
224 8 : CALL MPI_ALLREDUCE(this%dulo,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
225 8 : CALL zcopy(n,ctmp,1,this%dulo,1)
226 8 : CALL MPI_ALLREDUCE(this%ulod,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
227 8 : CALL zcopy(n,ctmp,1,this%ulod,1)
228 8 : DEALLOCATE(ctmp)
229 :
230 64 : n = SIZE(this%uloulop)
231 24 : ALLOCATE(ctmp(n))
232 8 : CALL MPI_ALLREDUCE(this%uloulop,ctmp,n,MPI_DOUBLE_COMPLEX,MPI_SUM,mpi_communicator,ierr)
233 8 : CALL zcopy(n,ctmp,1,this%uloulop,1)
234 8 : DEALLOCATE(ctmp)
235 : ENDIF
236 : ENDIF
237 : #endif
238 1018 : END SUBROUTINE collect_greensf
239 :
240 :
241 509 : FUNCTION occmtx_greensf_complete(this,gfinp,input,atoms,noco,nococonv,l_write,check,occError) Result(occmtx)
242 :
243 : !calculates the occupation of the greens function
244 : !The Greens-function should already be prepared on a energy contour ending at e_fermi
245 : !The occupation is calculated with:
246 : !
247 : ! n^sigma_mm' = -1/2pi int^Ef dz (G^+(z)^sigma_mm'-G^-(z)^sigma_mm')
248 : !
249 : ! If l_write is given the density matrix together with the spin up/down trace is written to the out files
250 :
251 : CLASS(t_greensf), INTENT(IN) :: this
252 : TYPE(t_gfinp), INTENT(IN) :: gfinp
253 : TYPE(t_input), INTENT(IN) :: input
254 : TYPE(t_atoms), INTENT(IN) :: atoms
255 : TYPE(t_noco), INTENT(IN) :: noco
256 : TYPE(t_nococonv), INTENT(IN) :: nococonv
257 : LOGICAL, OPTIONAL,INTENT(IN) :: l_write !write the occupation matrix to out file
258 : LOGICAL, OPTIONAL,INTENT(IN) :: check
259 : LOGICAL, OPTIONAL,INTENT(INOUT) :: occError
260 :
261 : COMPLEX, ALLOCATABLE :: occmtx(:,:,:)
262 :
263 : INTEGER :: nspins, ispin, m, mp
264 : LOGICAL :: all_occError, single_occError
265 : REAL :: nup, ndwn
266 : COMPLEX :: offd
267 : CHARACTER(len=2) :: l_type
268 : CHARACTER(len=8) :: l_form
269 : TYPE(t_contourInp) :: contourInp
270 :
271 509 : contourInp = gfinp%contour(this%elem%iContour)
272 :
273 509 : IF(ALLOCATED(this%gmmpMat)) nspins = SIZE(this%gmmpMat,4)
274 509 : IF(ALLOCATED(this%uu)) nspins = SIZE(this%uu,4)
275 :
276 59781 : ALLOCATE(occmtx(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const, nspins), source=cmplx_0)
277 :
278 509 : all_occError = .FALSE.
279 1531 : DO ispin = 1, nspins
280 : occmtx(:,:,ispin) = this%occmtx_greensf_spin(ispin,gfinp,input,atoms,noco,nococonv,&
281 58254 : check=check,occError=single_occError)
282 1531 : all_occError = all_occError.OR.single_occError
283 : ENDDO
284 509 : IF(PRESENT(occError)) occError = all_occError
285 :
286 : !Io-part (ATM this subroutine is only called from rank 0)
287 509 : IF(PRESENT(l_write)) THEN
288 509 : IF(l_write) THEN
289 : !Write to file
290 509 : WRITE (l_type,'(i2)') 2*(2*this%elem%l+1)
291 509 : l_form = '('//l_type//'f8.4)'
292 509 : IF(this%elem%isIntersite()) THEN
293 : 8998 FORMAT(/,"Occupation matrix obtained from the green's function for atom: ",I3," atomp: ",I3," atomDiff: ", 3f8.4," l: ",I3," lp: ",I3)
294 470 : WRITE(oUnit,8998) this%elem%atomType, this%elem%atomTypep, this%elem%atomDiff, this%elem%l, this%elem%lp
295 39 : ELSE IF (this%elem%isOffDiag()) THEN
296 : 8999 FORMAT(/,"Occupation matrix obtained from the green's function for atom: ",I3," l: ",I3," lp: ",I3)
297 6 : WRITE(oUnit,8999) this%elem%atomType, this%elem%l, this%elem%lp
298 : ELSE
299 : 9000 FORMAT(/,"Occupation matrix obtained from the green's function for atom: ",I3," l: ",I3)
300 33 : WRITE(oUnit,9000) this%elem%atomType, this%elem%l
301 : ENDIF
302 509 : WRITE(oUnit,"(A)") "In the |L,S> basis:"
303 1531 : DO ispin = 1, MERGE(3, input%jspins, gfinp%l_mperp)
304 1022 : WRITE(oUnit,'(A,I0)') "Spin: ", ispin
305 1531 : WRITE(oUnit,l_form) ((occmtx(m,mp,ispin),m=-this%elem%l, this%elem%l),mp=-this%elem%lp, this%elem%lp)
306 : ENDDO
307 :
308 509 : IF(this%elem%l.EQ.this%elem%lp) THEN
309 503 : nup = 0.0
310 3024 : DO m = -this%elem%l, this%elem%l
311 3024 : nup = nup + REAL(occmtx(m,m,1))
312 : ENDDO
313 503 : WRITE(oUnit,'(/,1x,A,I0,A,A,A,f8.4)') "l--> ",this%elem%l, " Contour(",TRIM(ADJUSTL(contourInp%label)),") Spin-Up trace: ", nup
314 :
315 503 : IF(input%jspins.EQ.2) THEN
316 503 : ndwn = 0.0
317 3024 : DO m = -this%elem%l, this%elem%l
318 3024 : ndwn = ndwn + REAL(occmtx(m,m,2))
319 : ENDDO
320 503 : WRITE(oUnit,'(1x,A,I0,A,A,A,f8.4)') "l--> ",this%elem%l, " Contour(",TRIM(ADJUSTL(contourInp%label)),") Spin-Down trace: ", ndwn
321 : ENDIF
322 :
323 503 : IF(gfinp%l_mperp) THEN
324 4 : offd = cmplx_0
325 24 : DO m = -this%elem%l, this%elem%l
326 24 : offd = offd + occmtx(m,m,3)
327 : ENDDO
328 4 : WRITE(oUnit,'(1x,A,I0,A,A,A,f8.4)') "l--> ",this%elem%l, " Contour(",TRIM(ADJUSTL(contourInp%label)),") Spin-Offd trace (x): ", REAL(offd)
329 4 : WRITE(oUnit,'(1x,A,I0,A,A,A,f8.4)') "l--> ",this%elem%l, " Contour(",TRIM(ADJUSTL(contourInp%label)),") Spin-Offd trace (y): ", AIMAG(offd)
330 : ENDIF
331 : ENDIF
332 : ENDIF
333 : ENDIF
334 :
335 509 : END FUNCTION occmtx_greensf_complete
336 :
337 2044 : FUNCTION occmtx_greensf_spin(this,spin,gfinp,input,atoms,noco,nococonv,check,occError) Result(occmtx)
338 :
339 : USE m_rotMMPmat
340 : USE m_types_mat
341 :
342 : !calculates the occupation of the greens function for a given spin
343 : !The Greens-function should already be prepared on a energy contour ending at e_fermi
344 : !The occupation is calculated with:
345 : !
346 : ! n^sigma_mm' = -1/2pi int^Ef dz (G^+(z)^sigma_mm'-G^-(z)^sigma_mm')
347 :
348 : CLASS(t_greensf), INTENT(IN) :: this
349 : INTEGER, INTENT(IN) :: spin
350 : TYPE(t_gfinp), INTENT(IN) :: gfinp
351 : TYPE(t_input), INTENT(IN) :: input
352 : TYPE(t_atoms), INTENT(IN) :: atoms
353 : TYPE(t_noco), INTENT(IN) :: noco
354 : TYPE(t_nococonv), INTENT(IN) :: nococonv
355 : LOGICAL, OPTIONAL,INTENT(IN) :: check
356 : LOGICAL, OPTIONAL,INTENT(INOUT) :: occError
357 :
358 : COMPLEX :: occmtx(-lmaxU_const:lmaxU_const,-lmaxU_const:lmaxU_const)
359 :
360 : INTEGER :: ind1,ind2,ipm,iz
361 : INTEGER :: l,lp,atomType,atomTypep,m,mp
362 : REAL :: tr
363 : COMPLEX :: weight
364 1022 : TYPE(t_mat) :: gmat
365 : CHARACTER(len=300) :: message
366 : TYPE(t_contourInp) :: contourInp
367 :
368 1022 : contourInp = gfinp%contour(this%elem%iContour)
369 1022 : IF(PRESENT(occError)) occError = .FALSE.
370 :
371 1022 : l = this%elem%l
372 1022 : lp = this%elem%lp
373 1022 : atomType = this%elem%atomType
374 1022 : atomTypep = this%elem%atomTypep
375 :
376 : !Check for Contours not reproducing occupations
377 1022 : IF(contourInp%shape.EQ.CONTOUR_SEMICIRCLE_CONST.AND.ABS(contourInp%et).GT.1e-12) &
378 0 : WRITE(oUnit,*) "Energy contour not ending at efermi: These are not the actual occupations"
379 1022 : IF(contourInp%shape.EQ.CONTOUR_DOS_CONST.AND..NOT.contourInp%l_dosfermi) &
380 2 : WRITE(oUnit,*) "Energy contour not weighted for occupations: These are not the actual occupations"
381 :
382 58254 : occmtx = cmplx_0
383 :
384 3066 : DO ipm = 1, 2
385 : !Integrate over the contour:
386 284944 : DO iz = 1, this%contour%nz
387 : !get the corresponding gf-matrix
388 282900 : weight = MERGE(this%contour%de(iz),conjg(this%contour%de(iz)),ipm.EQ.1)
389 282900 : CALL this%get(atoms,iz,ipm.EQ.2,spin,gmat)
390 282900 : ind1 = 0
391 1702516 : DO m = -l, l
392 1417572 : ind1 = ind1 + 1
393 1417572 : ind2 = 0
394 8822124 : DO mp = -lp,lp
395 7121652 : ind2 = ind2 + 1
396 : occmtx(m,mp) = occmtx(m,mp) + ImagUnit/tpi_const * (-1)**(ipm-1) * gmat%data_c(ind1,ind2) &
397 8539224 : * weight
398 : ENDDO
399 : ENDDO
400 : ENDDO
401 : !For the contour 3 (real Axis just shifted with sigma) we can add the tails on both ends
402 3066 : IF(contourInp%shape.EQ.CONTOUR_DOS_CONST.AND.contourInp%l_anacont) THEN
403 : !left tail
404 4 : weight = MERGE(this%contour%de(1),conjg(this%contour%de(1)),ipm.EQ.1)
405 4 : CALL this%get(atoms,1,ipm.EQ.2,spin,gmat)
406 4 : ind1 = 0
407 24 : DO m = -l, l
408 20 : ind1 = ind1 + 1
409 20 : ind2 = 0
410 124 : DO mp = -lp,lp
411 100 : ind2 = ind2 + 1
412 120 : occmtx(m,mp) = occmtx(m,mp) - 1/tpi_const * gmat%data_c(ind1,ind2) * weight
413 : ENDDO
414 : ENDDO
415 : !right tail
416 4 : weight = MERGE(this%contour%de(this%contour%nz),conjg(this%contour%de(this%contour%nz)),ipm.EQ.1)
417 4 : CALL this%get(atoms,this%contour%nz,ipm.EQ.2,spin,gmat)
418 4 : ind1 = 0
419 24 : DO m = -l, l
420 20 : ind1 = ind1 + 1
421 20 : ind2 = 0
422 124 : DO mp = -lp,lp
423 100 : ind2 = ind2 + 1
424 120 : occmtx(m,mp) = occmtx(m,mp) + 1/tpi_const * gmat%data_c(ind1,ind2) * weight
425 : ENDDO
426 : ENDDO
427 : ENDIF
428 : ENDDO
429 :
430 : !Rotate the occupation matrix into the global frame in real-space
431 1022 : IF(noco%l_noco) THEN
432 16 : IF (.NOT.this%elem%isIntersite()) THEN
433 912 : occmtx = rotMMPmat(occmtx,nococonv%alph(atomType),nococonv%beta(atomType),0.0,l)
434 : ENDIF
435 1006 : ELSE IF(noco%l_soc) THEN
436 456 : occmtx = rotMMPmat(occmtx,nococonv%phi,nococonv%theta,0.0,l)
437 : ENDIF
438 :
439 : !Sanity check are the occupations reasonable?
440 1022 : IF(PRESENT(check)) THEN
441 1022 : IF(check) THEN
442 56 : IF(spin<=input%jspins) THEN !Only the spin-diagonal parts
443 52 : tr = 0.0
444 340 : DO m = -l,l
445 288 : tr = tr + REAL(occmtx(m,m))/(3.0-input%jspins)
446 288 : IF(REAL(occmtx(m,m))/(3.0-input%jspins).GT. 1.05 .OR.&
447 340 : REAL(occmtx(m,m))/(3.0-input%jspins).LT.-0.01) THEN
448 :
449 0 : IF(PRESENT(occError)) THEN
450 0 : occError = .TRUE.
451 : ELSE
452 0 : WRITE(message,9100) spin,m,REAL(occmtx(m,m))
453 : 9100 FORMAT("Invalid element in mmpmat (spin ",I1,",m ",I2,"): ",f14.8)
454 0 : CALL juDFT_warn(TRIM(ADJUSTL(message)),calledby="occupationMatrix")
455 : ENDIF
456 : ENDIF
457 : ENDDO
458 52 : IF(tr.LT.-0.01.OR.tr.GT.2*l+1.1) THEN
459 0 : IF(PRESENT(occError)) THEN
460 0 : occError = .TRUE.
461 : ELSE
462 0 : WRITE(message,9110) spin,tr
463 : 9110 FORMAT("Invalid occupation for spin ",I1,": ",f14.8)
464 0 : CALL juDFT_warn(TRIM(ADJUSTL(message)),calledby="occupationMatrix")
465 : ENDIF
466 : ENDIF
467 : ENDIF
468 : ENDIF
469 : ENDIF
470 :
471 2044 : END FUNCTION occmtx_greensf_spin
472 :
473 : !----------------------------------------------------------------------------------
474 : ! Following this comment there are multiple definitions for functions
475 : ! to access the data in the greensFunction Type:
476 : ! get_gf -> Get the (m,mp) matrix of the spherically averaged GF
477 : ! at a certain energy point. If the correct scalar products
478 : ! are provided, the radial dependent GF can also be recombined here
479 : ! getRadial_gf -> Returns the radial and energy dependent GF for a certain spin
480 : ! and m,mp pair
481 : ! getRadialSpin_gf -> Returns the radial and energy dependent GF for a certain
482 : ! m,mp pair. Also returns the 2x2 spin matrix at that point
483 : ! set_gf -> Set the value of the (m,mp) Matrix at a
484 : ! certain energy point with an input matrix
485 : !----------------------------------------------------------------------------------
486 :
487 282908 : SUBROUTINE get_gf(this,atoms,iz,l_conjg,spin,gmat)
488 :
489 : USE m_types_mat
490 :
491 : !Returns the matrix belonging to energy point iz with l,lp,nType,nTypep
492 : !can also return the spherically averaged GF with the given scalar products
493 :
494 : CLASS(t_greensf), INTENT(IN) :: this
495 : TYPE(t_atoms), INTENT(IN) :: atoms
496 : INTEGER, INTENT(IN) :: iz
497 : LOGICAL, INTENT(IN) :: l_conjg
498 : INTEGER, INTENT(IN) :: spin
499 : TYPE(t_mat), INTENT(INOUT) :: gmat !Return matrix
500 :
501 : INTEGER matsize1,matsize2,ind1,ind2
502 : INTEGER m,mp,spin1,spin2,ipm,spin_start,spin_end,spin_ind,m_ind,mp_ind
503 : INTEGER l,lp,atomType,atomTypep,nspins,ilo,ilop,iLO_ind,iLOp_ind
504 :
505 : REAL :: uun,dun,udn,ddn
506 282908 : REAL :: uulon(atoms%nlod),dulon(atoms%nlod),ulodn(atoms%nlod),uloun(atoms%nlod)
507 282908 : REAL :: uloulopn(atoms%nlod,atoms%nlod)
508 :
509 282908 : IF(.NOT.this%l_calc) THEN
510 0 : CALL juDFT_error("The requested Green's Function element was not calculated", calledby="get_gf")
511 : ENDIF
512 :
513 282908 : l = this%elem%l
514 282908 : lp = this%elem%lp
515 282908 : atomType = this%elem%atomType
516 282908 : atomTypep = this%elem%atomTypep
517 :
518 282908 : IF(ALLOCATED(this%gmmpMat)) THEN
519 279836 : nspins = SIZE(this%gmmpMat,4)
520 : ELSE
521 3072 : nspins = SIZE(this%uu,4)
522 : ENDIF
523 :
524 282908 : IF(spin.GT.3 .OR. spin.LT.1) THEN
525 0 : CALL juDFT_error("Invalid argument for spin",calledby="get_gf")
526 : ENDIF
527 :
528 282908 : matsize1 = 2*l+1
529 282908 : matsize2 = 2*lp+1
530 :
531 282908 : IF(.NOT.ALLOCATED(gmat%data_c)) THEN
532 1022 : CALL gmat%init(.FALSE.,matsize1,matsize2)
533 281886 : ELSE IF(matsize1.NE.gmat%matsize1.OR.matsize2.NE.gmat%matsize2) THEN
534 0 : CALL juDFT_error("Mismatch in matsizes", calledby="get_gf")
535 : ENDIF
536 :
537 282908 : ipm = MERGE(2,1,l_conjg)
538 :
539 8822372 : gmat%data_c = cmplx_0
540 :
541 282908 : IF(spin < 3) THEN
542 281884 : spin1 = spin
543 281884 : spin2 = spin
544 : ELSE
545 : spin1 = 2
546 : spin2 = 1
547 : ENDIF
548 :
549 : !Find the correct spin index in gmmpMat arrays
550 282908 : spin_ind = MERGE(1,spin,nspins.EQ.1)
551 :
552 282908 : IF(.NOT.this%l_sphavg) THEN
553 3072 : IF(.NOT.ALLOCATED(this%scalarProducts%uun)) THEN
554 0 : CALL juDFT_error('Scalar products not available')
555 : ENDIF
556 3072 : uun = this%scalarProducts%uun(spin1,spin2)
557 3072 : dun = this%scalarProducts%dun(spin1,spin2)
558 3072 : udn = this%scalarProducts%udn(spin1,spin2)
559 3072 : ddn = this%scalarProducts%ddn(spin1,spin2)
560 :
561 3072 : IF(ALLOCATED(this%uulo)) THEN
562 8704 : uulon(:) = this%scalarProducts%uulon(:,spin1,spin2)
563 8704 : uloun(:) = this%scalarProducts%uloun(:,spin1,spin2)
564 8704 : dulon(:) = this%scalarProducts%dulon(:,spin1,spin2)
565 8704 : ulodn(:) = this%scalarProducts%ulodn(:,spin1,spin2)
566 :
567 31744 : uloulopn(:,:) = this%scalarProducts%uloulopn(:,:,spin1,spin2)
568 : ENDIF
569 : ENDIF
570 :
571 282908 : ind1 = 0
572 1700520 : DO m = -l,l
573 1417612 : ind1 = ind1 + 1
574 1417612 : ind2 = 0
575 8822372 : DO mp = -lp,lp
576 7121852 : ind2 = ind2 + 1
577 :
578 : !-------------------------------------------------------------------
579 : ! Check wether we need to do some operation on the indices m and mp
580 : !-------------------------------------------------------------------
581 7121852 : IF(spin.EQ.2 .AND. nspins.EQ.1) THEN
582 : !For a non-spin-polarized calculation we might still want the full
583 : !matrix. Then we need to reverse the order (SOC prop m*s_z)
584 0 : m_ind = -m
585 0 : mp_ind = -mp
586 : ELSE
587 : !Do nothing
588 : m_ind = m
589 : mp_ind = mp
590 : ENDIF
591 : !-------------------
592 : ! Fetch the values
593 : !-------------------
594 8539464 : IF(.NOT.this%l_sphavg) THEN
595 : gmat%data_c(ind1,ind2) = this%uu(iz,m_ind,mp_ind,spin_ind,ipm) * uun &
596 : + this%dd(iz,m_ind,mp_ind,spin_ind,ipm) * ddn &
597 : + this%du(iz,m_ind,mp_ind,spin_ind,ipm) * dun &
598 60416 : + this%ud(iz,m_ind,mp_ind,spin_ind,ipm) * udn
599 60416 : IF(ALLOCATED(this%uulo)) THEN
600 34816 : iLO_ind = 0
601 152064 : DO ilo = 1, atoms%nlo(atomType)
602 117248 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
603 39424 : iLO_ind = iLO_ind + 1
604 : gmat%data_c(ind1,ind2) = gmat%data_c(ind1,ind2) &
605 : + this%uulo(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * uulon(ilo) &
606 152064 : + this%dulo(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * dulon(ilo)
607 : ENDDO
608 34816 : iLO_ind = 0
609 152064 : DO ilo = 1, atoms%nlo(atomTypep)
610 117248 : IF(atoms%llo(ilo,atomTypep).NE.lp) CYCLE
611 39424 : iLO_ind = iLO_ind + 1
612 : gmat%data_c(ind1,ind2) = gmat%data_c(ind1,ind2) &
613 : + this%ulou(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * uloun(ilo) &
614 152064 : + this%dulo(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ulodn(ilo)
615 : ENDDO
616 : iLO_ind = 0
617 152064 : DO ilo = 1,atoms%nlo(atomType)
618 117248 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
619 39424 : iLO_ind = iLO_ind + 1
620 39424 : iLOp_ind = 0
621 209920 : DO ilop = 1, atoms%nlo(atomTypep)
622 135680 : IF(atoms%llo(ilop,atomTypep).NE.lp) CYCLE
623 48640 : iLOp_ind = iLOp_ind + 1
624 : gmat%data_c(ind1,ind2) = gmat%data_c(ind1,ind2) &
625 : + this%uloulop(iz,m_ind,mp_ind,iLO_ind,iLOp_ind,spin_ind,ipm) &
626 252928 : * uloulopn(ilo,ilop)
627 : ENDDO
628 : ENDDO
629 : ENDIF
630 : ELSE
631 7061436 : gmat%data_c(ind1,ind2) = this%gmmpMat(iz,m_ind,mp_ind,spin_ind,ipm)
632 : ENDIF
633 :
634 : ENDDO!mp
635 : ENDDO!m
636 :
637 282908 : END SUBROUTINE get_gf
638 :
639 0 : SUBROUTINE getFullMatrix_gf(this,atoms,iz,l_conjg,gmat)
640 :
641 : USE m_types_mat
642 :
643 : !Return the full matrix with all spin blocks for the given energy point
644 :
645 : CLASS(t_greensf), INTENT(IN) :: this
646 : TYPE(t_atoms), INTENT(IN) :: atoms
647 : INTEGER, INTENT(IN) :: iz
648 : LOGICAL, INTENT(IN) :: l_conjg
649 : TYPE(t_mat), INTENT(INOUT) :: gmat !Return matrix
650 :
651 : INTEGER :: matsize1, matsize2, nspins, ispin
652 :
653 0 : TYPE(t_mat) :: gmat_spin
654 :
655 0 : matsize1 = 2*this%elem%l+1
656 0 : matsize2 = 2*this%elem%l+1
657 :
658 0 : IF(.NOT.ALLOCATED(gmat%data_c)) THEN
659 0 : CALL gmat%init(.FALSE.,2*matsize1,2*matsize2)
660 0 : ELSE IF(2*matsize1.NE.gmat%matsize1.OR.2*matsize2.NE.gmat%matsize2) THEN
661 0 : CALL juDFT_error("Mismatch in matsizes", calledby="getFullMatrix_gf")
662 : ENDIF
663 :
664 0 : IF(ALLOCATED(this%gmmpMat)) THEN
665 0 : nspins = SIZE(this%gmmpMat,4)
666 : ELSE
667 0 : nspins = SIZE(this%uu,4)
668 : ENDIF
669 :
670 0 : DO ispin = 1, MAX(nspins,2)
671 0 : CALL this%get(atoms,iz,l_conjg,ispin,gmat_spin)
672 :
673 0 : IF(ispin<3) THEN
674 0 : gmat%data_c((ispin-1)*matsize1+1:ispin*matsize1,(ispin-1)*matsize2+1:ispin*matsize2) = gmat_spin%data_c
675 : ELSE
676 0 : gmat%data_c(matsize1+1:2*matsize1,1:matsize2) = gmat_spin%data_c
677 0 : gmat%data_c(1:matsize1,matsize2+1:2*matsize2) = conjg(transpose(gmat_spin%data_c))
678 : ENDIF
679 : ENDDO
680 :
681 0 : IF(nspins==1) gmat%data_c = gmat%data_c * 0.5
682 :
683 0 : END SUBROUTINE getFullMatrix_gf
684 :
685 0 : SUBROUTINE getRadial_gf(this,atoms,m,mp,l_conjg,spin,f,g,flo,gmat)
686 :
687 : !Returns the green's function on the radial and energy mesh
688 : !for a certain m,mp,spin combination. Attention: The correct radial functions have to be provided
689 :
690 : CLASS(t_greensf), INTENT(IN) :: this
691 : TYPE(t_atoms), INTENT(IN) :: atoms
692 : INTEGER, INTENT(IN) :: m,mp
693 : LOGICAL, INTENT(IN) :: l_conjg
694 : INTEGER, INTENT(IN) :: spin
695 : REAL , INTENT(IN) :: f(:,:,0:,:,:)
696 : REAL , INTENT(IN) :: g(:,:,0:,:,:)
697 : REAL , INTENT(IN) :: flo(:,:,:,:,:)
698 : COMPLEX, ALLOCATABLE,INTENT(INOUT) :: gmat(:,:) !Return matrix
699 :
700 : INTEGER spin1,spin2,ipm,spin_ind,m_ind,mp_ind,ilo,ilop,iLO_ind,iLOp_ind
701 : INTEGER l,lp,atomType,atomTypep,nspins,iz
702 :
703 0 : IF(.NOT.this%l_calc) THEN
704 0 : CALL juDFT_error("The requested Green's Function element was not calculated", calledby="get_gf")
705 : ENDIF
706 :
707 0 : l = this%elem%l
708 0 : lp = this%elem%lp
709 0 : atomType = this%elem%atomType
710 0 : atomTypep = this%elem%atomTypep
711 :
712 0 : IF(ALLOCATED(this%gmmpMat)) THEN
713 0 : CALL juDFT_error("Green's function not calculated for radial dependence", calledby="get_gf")
714 : ENDIF
715 :
716 0 : nspins = SIZE(this%uu,4)
717 :
718 0 : IF(spin.GT.3 .OR. spin.LT.1) THEN
719 0 : CALL juDFT_error("Invalid argument for spin",calledby="get_gf")
720 : ENDIF
721 :
722 0 : ipm = MERGE(2,1,l_conjg)
723 :
724 0 : IF(.NOT.ALLOCATED(gmat)) ALLOCATE(gmat(SIZE(f,1),this%contour%nz),source=cmplx_0)
725 0 : gmat = cmplx_0
726 :
727 0 : IF(spin < 3) THEN
728 0 : spin1 = spin
729 0 : spin2 = spin
730 : ELSE
731 : spin1 = 2
732 : spin2 = 1
733 : ENDIF
734 : !Find the correct spin index in gmmpMat arrays
735 0 : spin_ind = MERGE(1,spin,nspins.EQ.1)
736 :
737 : !-------------------------------------------------------------------
738 : ! Check wether we need to do some operation on the indices m and mp
739 : !-------------------------------------------------------------------
740 0 : IF(spin.EQ.2 .AND. nspins.EQ.1) THEN
741 : !For a non-spin-polarized calculation we might still want the full
742 : !matrix. Then we need to reverse the order (SOC prop m*s_z)
743 0 : m_ind = -m
744 0 : mp_ind = -mp
745 : ELSE
746 : !Do nothing
747 0 : m_ind = m
748 0 : mp_ind = mp
749 : ENDIF
750 : !-------------------
751 : ! Fetch the values
752 : !-------------------
753 0 : DO iz = 1, this%contour%nz
754 : gmat(:,iz) = this%uu(iz,m_ind,mp_ind,spin_ind,ipm) * ( f(:,1,l,spin2,atomType) * f(:,1,lp,spin1,atomTypep) &
755 : +f(:,2,l,spin2,atomType) * f(:,2,lp,spin1,atomTypep))&
756 : + this%dd(iz,m_ind,mp_ind,spin_ind,ipm) * ( g(:,1,l,spin2,atomType) * g(:,1,lp,spin1,atomTypep) &
757 : +g(:,2,l,spin2,atomType) * g(:,2,lp,spin1,atomTypep))&
758 : + this%ud(iz,m_ind,mp_ind,spin_ind,ipm) * ( g(:,1,l,spin2,atomType) * f(:,1,lp,spin1,atomTypep) &
759 : +g(:,2,l,spin2,atomType) * f(:,2,lp,spin1,atomTypep))&
760 : + this%du(iz,m_ind,mp_ind,spin_ind,ipm) * ( f(:,1,l,spin2,atomType) * g(:,1,lp,spin1,atomTypep) &
761 0 : +f(:,2,l,spin2,atomType) * g(:,2,lp,spin1,atomTypep))
762 :
763 0 : IF(ALLOCATED(this%uulo)) THEN
764 0 : iLO_ind = 0
765 0 : DO ilo = 1, atoms%nlo(atomType)
766 0 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
767 0 : iLO_ind = iLO_ind + 1
768 : gmat(:,iz) = gmat(:,iz) &
769 : + this%uulo(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ( f(:,1,lp,spin1,atomTypep) *flo(:,1,ilo,spin2,atomType) &
770 : +f(:,2,lp,spin1,atomTypep) *flo(:,2,ilo,spin2,atomType))&
771 : + this%dulo(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ( g(:,1,lp,spin1,atomTypep) *flo(:,1,ilo,spin2,atomType) &
772 0 : +g(:,2,lp,spin1,atomTypep) *flo(:,2,ilo,spin2,atomType))
773 : ENDDO
774 0 : iLO_ind = 0
775 0 : DO ilo = 1, atoms%nlo(atomTypep)
776 0 : IF(atoms%llo(ilo,atomTypep).NE.lp) CYCLE
777 0 : iLO_ind = iLO_ind + 1
778 : gmat(:,iz) = gmat(:,iz) &
779 : + this%ulou(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ( flo(:,1,ilo,spin1,atomTypep)*f(:,1,l,spin2,atomType) &
780 : +flo(:,2,ilo,spin1,atomTypep)*f(:,2,l,spin2,atomType))&
781 : + this%ulod(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ( flo(:,1,ilo,spin1,atomTypep)*g(:,1,l,spin2,atomType) &
782 0 : +flo(:,2,ilo,spin1,atomTypep)*g(:,2,l,spin2,atomType))
783 : ENDDO
784 : iLO_ind = 0
785 0 : DO ilo = 1, atoms%nlo(atomType)
786 0 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
787 0 : iLOp_ind = 0
788 0 : iLO_ind = iLO_ind + 1
789 0 : DO ilop = 1, atoms%nlo(atomTypep)
790 0 : IF(atoms%llo(ilop,atomTypep).NE.lp) CYCLE
791 0 : iLOp_ind = iLOp_ind + 1
792 : gmat(:,iz) = gmat(:,iz) &
793 : + this%uloulop(iz,m_ind,mp_ind,iLO_ind,iLOp_ind,spin_ind,ipm) *( flo(:,1,ilo,spin2,atomType)*flo(:,1,ilop,spin1,atomTypep) &
794 0 : +flo(:,2,ilo,spin2,atomType)*flo(:,2,ilop,spin1,atomTypep))
795 : ENDDO
796 : ENDDO
797 : ENDIF
798 : ENDDO
799 :
800 0 : END SUBROUTINE getRadial_gf
801 :
802 0 : SUBROUTINE getRadialRadial_gf(this,atoms,iz,m,mp,l_conjg,spin,f,g,flo,gmat)
803 :
804 : !Returns the green's function on the radial and energy mesh (r/=r')
805 : !for a certain m,mp,spin combination. Attention: The correct radial functions have to be provided
806 :
807 : CLASS(t_greensf), INTENT(IN) :: this
808 : TYPE(t_atoms), INTENT(IN) :: atoms
809 : INTEGER, INTENT(IN) :: iz
810 : INTEGER, INTENT(IN) :: m,mp
811 : LOGICAL, INTENT(IN) :: l_conjg
812 : INTEGER, INTENT(IN) :: spin
813 : REAL , INTENT(IN) :: f(:,:,0:,:,:)
814 : REAL , INTENT(IN) :: g(:,:,0:,:,:)
815 : REAL , INTENT(IN) :: flo(:,:,:,:,:)
816 : COMPLEX, ALLOCATABLE,INTENT(INOUT) :: gmat(:,:) !Return matrix
817 :
818 : INTEGER spin1,spin2,ipm,spin_ind,m_ind,mp_ind,ilo,ilop,iLO_ind,iLOp_ind
819 : INTEGER l,lp,atomType,atomTypep,nspins,jr,jrp
820 :
821 0 : IF(.NOT.this%l_calc) THEN
822 0 : CALL juDFT_error("The requested Green's Function element was not calculated", calledby="get_gf")
823 : ENDIF
824 :
825 0 : l = this%elem%l
826 0 : lp = this%elem%lp
827 0 : atomType = this%elem%atomType
828 0 : atomTypep = this%elem%atomTypep
829 :
830 0 : IF(ALLOCATED(this%gmmpMat)) THEN
831 0 : CALL juDFT_error("Green's function not calculated for radial dependence", calledby="get_gf")
832 : ENDIF
833 :
834 0 : nspins = SIZE(this%uu,4)
835 :
836 0 : IF(spin.GT.3 .OR. spin.LT.1) THEN
837 0 : CALL juDFT_error("Invalid argument for spin",calledby="get_gf")
838 : ENDIF
839 :
840 0 : ipm = MERGE(2,1,l_conjg)
841 :
842 0 : IF(.NOT.ALLOCATED(gmat)) ALLOCATE(gmat(SIZE(f,1),SIZE(f,1)),source=cmplx_0)
843 :
844 0 : IF(spin < 3) THEN
845 0 : spin1 = spin
846 0 : spin2 = spin
847 : ELSE
848 : spin1 = 2
849 : spin2 = 1
850 : ENDIF
851 : !Find the correct spin index in gmmpMat arrays
852 0 : spin_ind = MERGE(1,spin,nspins.EQ.1)
853 :
854 : !-------------------------------------------------------------------
855 : ! Check wether we need to do some operation on the indices m and mp
856 : !-------------------------------------------------------------------
857 0 : IF(spin.EQ.2 .AND. nspins.EQ.1) THEN
858 : !For a non-spin-polarized calculation we might still want the full
859 : !matrix. Then we need to reverse the order (SOC prop m*s_z)
860 0 : m_ind = -m
861 0 : mp_ind = -mp
862 : ELSE
863 : !Do nothing
864 0 : m_ind = m
865 0 : mp_ind = mp
866 : ENDIF
867 : !-------------------
868 : ! Fetch the values
869 : !-------------------
870 : !$OMP parallel do default(none) &
871 : !$OMP shared(this,atoms,f,g,flo,gmat,m_ind,mp_ind,spin_ind,ipm,atomType,atomTypep,spin1,spin2,l,lp,iz) &
872 0 : !$OMP private(jr,iLO_ind,iLOp_ind,ilo,ilop)
873 : DO jr = 1, atoms%jri(atomType)
874 : gmat(:,jr) = this%uu(iz,m_ind,mp_ind,spin_ind,ipm) * ( f(jr,1,l,spin2,atomType) * f(:,1,lp,spin1,atomTypep) &
875 : +f(jr,2,l,spin2,atomType) * f(:,2,lp,spin1,atomTypep))&
876 : + this%dd(iz,m_ind,mp_ind,spin_ind,ipm) * ( g(jr,1,l,spin2,atomType) * g(:,1,lp,spin1,atomTypep) &
877 : +g(jr,2,l,spin2,atomType) * g(:,2,lp,spin1,atomTypep))&
878 : + this%ud(iz,m_ind,mp_ind,spin_ind,ipm) * ( g(jr,1,l,spin2,atomType) * f(:,1,lp,spin1,atomTypep) &
879 : +g(jr,2,l,spin2,atomType) * f(:,2,lp,spin1,atomTypep))&
880 : + this%du(iz,m_ind,mp_ind,spin_ind,ipm) * ( f(jr,1,l,spin2,atomType) * g(:,1,lp,spin1,atomTypep) &
881 : +f(jr,2,l,spin2,atomType) * g(:,2,lp,spin1,atomTypep))
882 :
883 : IF(ALLOCATED(this%uulo)) THEN
884 : iLO_ind = 0
885 : DO ilo = 1, atoms%nlo(atomType)
886 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
887 : iLO_ind = iLO_ind + 1
888 : gmat(:,jr) = gmat(:,jr) &
889 : + this%uulo(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ( f(:,1,lp,spin1,atomTypep) *flo(jr,1,ilo,spin2,atomType) &
890 : +f(:,2,lp,spin1,atomTypep) *flo(jr,2,ilo,spin2,atomType))&
891 : + this%dulo(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ( g(:,1,lp,spin1,atomTypep) *flo(jr,1,ilo,spin2,atomType) &
892 : +g(:,2,lp,spin1,atomTypep) *flo(jr,2,ilo,spin2,atomType))
893 : ENDDO
894 : iLO_ind = 0
895 : DO ilo = 1, atoms%nlo(atomTypep)
896 : IF(atoms%llo(ilo,atomTypep).NE.lp) CYCLE
897 : iLO_ind = iLO_ind + 1
898 : gmat(:,jr) = gmat(:,jr) &
899 : + this%ulou(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ( flo(:,1,ilo,spin1,atomTypep)*f(jr,1,l,spin2,atomType) &
900 : +flo(:,2,ilo,spin1,atomTypep)*f(jr,2,l,spin2,atomType))&
901 : + this%ulod(iz,m_ind,mp_ind,iLO_ind,spin_ind,ipm) * ( flo(:,1,ilo,spin1,atomTypep)*g(jr,1,l,spin2,atomType) &
902 : +flo(:,2,ilo,spin1,atomTypep)*g(jr,2,l,spin2,atomType))
903 : ENDDO
904 : iLO_ind = 0
905 : DO ilo = 1, atoms%nlo(atomType)
906 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
907 : iLOp_ind = 0
908 : iLO_ind = iLO_ind + 1
909 : DO ilop = 1, atoms%nlo(atomTypep)
910 : IF(atoms%llo(ilop,atomType).NE.lp) CYCLE
911 : iLOp_ind = iLOp_ind + 1
912 : gmat(:,jr) = gmat(:,jr) &
913 : + this%uloulop(iz,m_ind,mp_ind,iLO_ind,iLOp_ind,spin_ind,ipm) *( flo(jr,1,ilo,spin2,atomType)*flo(:,1,ilop,spin1,atomTypep) &
914 : +flo(jr,2,ilo,spin2,atomType)*flo(:,2,ilop,spin1,atomTypep))
915 : ENDDO
916 : ENDDO
917 : ENDIF
918 : !Get the right normalization
919 : gmat(:,jr) = gmat(:,jr) * atoms%rmsh(jr,atomType) * atoms%rmsh(:atoms%jri(atomTypep),atomTypep)
920 : ENDDO
921 : !$OMP end parallel do
922 :
923 0 : END SUBROUTINE getRadialRadial_gf
924 :
925 0 : SUBROUTINE set_gf(this,iz,l_conjg,gmat,spin)
926 :
927 : USE m_types_mat
928 :
929 : !Sets the spherically averaged greens function matrix belonging to energy point iz with l,lp,nType,nTypep
930 : !equal to gmat
931 :
932 : CLASS(t_greensf), INTENT(INOUT) :: this
933 : INTEGER, INTENT(IN) :: iz
934 : LOGICAL, INTENT(IN) :: l_conjg
935 : TYPE(t_mat), INTENT(IN) :: gmat
936 : INTEGER, OPTIONAL, INTENT(IN) :: spin
937 :
938 : INTEGER matsize1,matsize2,ind1,ind2,ind1_start,ind2_start
939 : INTEGER l,lp,atomType,atomTypep,m,mp,spin1,spin2,ipm,ispin,spin_start,spin_end
940 : LOGICAL l_full
941 :
942 :
943 0 : this%l_calc = .TRUE. !If its set it counts as calculated
944 :
945 0 : l = this%elem%l
946 0 : lp = this%elem%lp
947 0 : atomType = this%elem%atomType
948 0 : atomTypep = this%elem%atomTypep
949 :
950 :
951 0 : IF(ALLOCATED(this%uu)) THEN
952 0 : CALL juDFT_error("Can only set spherically averaged Green's Function", calledby="set_gf")
953 : ENDIF
954 :
955 0 : IF(PRESENT(spin)) THEN
956 0 : IF(spin.GT.3 .OR. spin.LT.1) THEN
957 0 : CALL juDFT_error("Invalid argument for spin",calledby="set_gf")
958 : ENDIF
959 : ENDIF
960 :
961 0 : l_full = .NOT.PRESENT(spin)
962 : !Determine matsize for the result gmat (if spin is given only return this digonal element)
963 0 : matsize1 = (2*l+1) * MERGE(2,1,l_full)
964 0 : matsize2 = (2*lp+1) * MERGE(2,1,l_full)
965 :
966 : !Check the expected matsizes against the actual
967 0 : IF(matsize1.NE.gmat%matsize1.OR.matsize2.NE.gmat%matsize2) THEN
968 0 : CALL juDFT_error("Mismatch in matsizes", calledby="set_gf")
969 : ENDIF
970 :
971 0 : ipm = MERGE(2,1,l_conjg)
972 :
973 0 : IF(l_full) THEN
974 0 : spin_start = 1
975 0 : spin_end = SIZE(this%gmmpMat,4)
976 : ELSE
977 0 : spin_start = spin
978 0 : spin_end = spin
979 : ENDIF
980 :
981 0 : DO ispin = spin_start, spin_end
982 : !Find the right quadrant in gmat according to the spin index
983 0 : IF(ispin.EQ.2 .AND.SIZE(this%gmmpMat,4).EQ.1) CYCLE
984 0 : IF(l_full) THEN
985 0 : IF(ispin < 3) THEN
986 0 : spin1 = ispin
987 0 : spin2 = ispin
988 : ELSE
989 : spin1 = 2
990 : spin2 = 1
991 : ENDIF
992 0 : ind1_start = (spin1-1)*(2*l+1)
993 0 : ind2_start = (spin2-1)*(2*lp+1)
994 : ELSE
995 : ind1_start = 0
996 : ind2_start = 0
997 : ENDIF
998 0 : ind1 = ind1_start
999 0 : DO m = -l,l
1000 0 : ind1 = ind1 + 1
1001 0 : ind2 = ind2_start
1002 0 : DO mp = -lp,lp
1003 0 : ind2 = ind2 + 1
1004 0 : this%gmmpMat(iz,m,mp,ispin,ipm) = gmat%data_c(ind1,ind2)
1005 0 : IF(l_full) this%gmmpMat(iz,m,mp,ispin,ipm) = this%gmmpMat(iz,m,mp,ispin,ipm) &
1006 0 : * MERGE(2.0,1.0,SIZE(this%gmmpMat,4).EQ.1)
1007 : ENDDO
1008 : ENDDO
1009 : ENDDO
1010 :
1011 0 : END SUBROUTINE set_gf
1012 :
1013 377 : SUBROUTINE set_gfdata(this,repr_gf)
1014 :
1015 : !Sets all arrays equal to the data from another greensfunction
1016 :
1017 : CLASS(t_greensf), INTENT(INOUT) :: this
1018 : TYPE(t_greensf), INTENT(IN) :: repr_gf
1019 :
1020 : !TODO check array sizes before
1021 :
1022 9545640 : IF(ALLOCATED(repr_gf%gmmpMat)) this%gmmpMat = repr_gf%gmmpMat
1023 377 : IF(ALLOCATED(repr_gf%gmmpMat_k)) this%gmmpMat_k = repr_gf%gmmpMat_k
1024 377 : IF(ALLOCATED(repr_gf%uu)) this%uu = repr_gf%uu
1025 377 : IF(ALLOCATED(repr_gf%ud)) this%ud = repr_gf%ud
1026 377 : IF(ALLOCATED(repr_gf%du)) this%du = repr_gf%du
1027 377 : IF(ALLOCATED(repr_gf%dd)) this%dd = repr_gf%dd
1028 377 : IF(ALLOCATED(repr_gf%uulo)) this%uulo = repr_gf%uulo
1029 377 : IF(ALLOCATED(repr_gf%ulou)) this%ulou = repr_gf%ulou
1030 377 : IF(ALLOCATED(repr_gf%ulod)) this%ulod = repr_gf%ulod
1031 377 : IF(ALLOCATED(repr_gf%dulo)) this%dulo = repr_gf%dulo
1032 377 : IF(ALLOCATED(repr_gf%uloulop)) this%uloulop = repr_gf%uloulop
1033 :
1034 377 : END SUBROUTINE set_gfdata
1035 :
1036 377 : SUBROUTINE rotate_gf(this,sym,atoms)
1037 :
1038 : USE m_rotMMPmat
1039 :
1040 : !Applies the given symmetry operation to the greens function
1041 : CLASS(t_greensf), INTENT(INOUT) :: this
1042 : TYPE(t_sym), INTENT(IN) :: sym
1043 : TYPE(t_atoms), INTENT(IN) :: atoms
1044 :
1045 : INTEGER :: l,lp,iop,atomType,atomTypep,ikpt
1046 : INTEGER :: ipm,ispin,iz,ilo,ilop,iLO_ind,iLOp_ind
1047 :
1048 377 : IF(this%elem%representative_elem<0) RETURN !Nothing to be done
1049 :
1050 377 : CALL timestart("Green's Function: Rotate (Symmetry)")
1051 377 : l = this%elem%l
1052 377 : lp = this%elem%lp
1053 377 : atomType = this%elem%atomType
1054 377 : atomTypep = this%elem%atomTypep
1055 377 : iop = this%elem%representative_op
1056 :
1057 1131 : DO ipm = 1, 2
1058 97643 : DO iz = 1, this%contour%nz
1059 97266 : IF(ALLOCATED(this%gmmpMat)) THEN
1060 11098880 : this%gmmpMat(iz,:,:,:,ipm) = rotMMPmat(this%gmmpMat(iz,:,:,:,ipm),sym,iop,l,lp=lp)
1061 0 : ELSE IF(ALLOCATED(this%uu)) THEN
1062 0 : this%uu(iz,:,:,:,ipm) = rotMMPmat(this%uu(iz,:,:,:,ipm),sym,iop,l,lp=lp)
1063 0 : this%dd(iz,:,:,:,ipm) = rotMMPmat(this%dd(iz,:,:,:,ipm),sym,iop,l,lp=lp)
1064 0 : this%ud(iz,:,:,:,ipm) = rotMMPmat(this%ud(iz,:,:,:,ipm),sym,iop,l,lp=lp)
1065 0 : this%du(iz,:,:,:,ipm) = rotMMPmat(this%du(iz,:,:,:,ipm),sym,iop,l,lp=lp)
1066 :
1067 0 : IF(ALLOCATED(this%uulo)) THEN
1068 0 : iLO_ind = 0
1069 0 : DO ilo = 1, atoms%nlo(atomType)
1070 0 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
1071 0 : iLO_ind = iLO_ind + 1
1072 0 : this%uulo(iz,:,:,iLO_ind,:,ipm) = rotMMPmat(this%uulo(iz,:,:,iLO_ind,:,ipm),sym,iop,l,lp=lp)
1073 0 : this%dulo(iz,:,:,iLO_ind,:,ipm) = rotMMPmat(this%dulo(iz,:,:,iLO_ind,:,ipm),sym,iop,l,lp=lp)
1074 : ENDDO
1075 0 : iLO_ind = 0
1076 0 : DO ilo = 1, atoms%nlo(atomTypep)
1077 0 : IF(atoms%llo(ilo,atomTypep).NE.lp) CYCLE
1078 0 : iLO_ind = iLO_ind + 1
1079 0 : this%ulou(iz,:,:,iLO_ind,:,ipm) = rotMMPmat(this%ulou(iz,:,:,iLO_ind,:,ipm),sym,iop,l,lp=lp)
1080 0 : this%ulod(iz,:,:,iLO_ind,:,ipm) = rotMMPmat(this%ulod(iz,:,:,iLO_ind,:,ipm),sym,iop,l,lp=lp)
1081 : ENDDO
1082 0 : iLO_ind = 0
1083 0 : DO ilo = 1, atoms%nlo(atomType)
1084 0 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
1085 0 : iLOp_ind = 0
1086 0 : iLO_ind = iLO_ind + 1
1087 0 : DO ilop = 1, atoms%nlo(atomTypep)
1088 0 : IF(atoms%llo(ilop,atomType).NE.lp) CYCLE
1089 0 : iLOp_ind = iLOp_ind + 1
1090 0 : this%uloulop(iz,:,:,iLO_ind,iLOp_ind,:,ipm) = rotMMPmat(this%uloulop(iz,:,:,iLO_ind,iLOp_ind,:,ipm),sym,iop,l,lp=lp)
1091 : ENDDO
1092 : ENDDO
1093 : ENDIF
1094 0 : ELSE IF(ALLOCATED(this%gmmpMat_k)) THEN
1095 0 : DO ikpt = 1, SIZE(this%gmmpMat_k,6)
1096 0 : this%gmmpMat_k(iz,:,:,:,ipm,ikpt) = rotMMPmat(this%gmmpMat_k(iz,:,:,:,ipm,ikpt),sym,iop,l,lp=lp)
1097 : ENDDO
1098 : ENDIF
1099 : ENDDO
1100 : ENDDO
1101 377 : CALL timestop("Green's Function: Rotate (Symmetry)")
1102 : END SUBROUTINE rotate_gf
1103 :
1104 0 : SUBROUTINE rotate_euler_angles_gf(this,atoms,alpha,beta,gamma,spin_rotation,real_space_rotation)
1105 :
1106 : USE m_rotMMPmat
1107 :
1108 : !Applies the given symmetry operation to the greens function
1109 : CLASS(t_greensf), INTENT(INOUT) :: this
1110 : TYPE(t_atoms), INTENT(IN) :: atoms
1111 : REAL, INTENT(IN) :: alpha
1112 : REAL, INTENT(IN) :: beta
1113 : REAL, INTENT(IN) :: gamma
1114 : LOGICAL, OPTIONAL, INTENT(IN) :: spin_rotation
1115 : LOGICAL, OPTIONAL, INTENT(IN) :: real_space_rotation
1116 :
1117 : INTEGER :: l,lp,atomType,atomTypep,ikpt
1118 : INTEGER :: ipm,ispin,iz,ilo,ilop,iLO_ind,iLOp_ind
1119 :
1120 0 : IF(ABS(alpha).LT.1e-12.AND.ABS(beta).LT.1e-12.AND.ABS(gamma).LT.1e-12) RETURN !Nothing to be done
1121 :
1122 0 : CALL timestart("Green's Function: Rotate (Angles)")
1123 0 : l = this%elem%l
1124 0 : lp = this%elem%lp
1125 0 : atomType = this%elem%atomType
1126 0 : atomTypep = this%elem%atomTypep
1127 :
1128 0 : DO ipm = 1, 2
1129 0 : DO iz = 1, this%contour%nz
1130 0 : IF(ALLOCATED(this%gmmpMat)) THEN
1131 : this%gmmpMat(iz,:,:,:,ipm) = rotMMPmat(this%gmmpMat(iz,:,:,:,ipm),alpha,beta,gamma,l,lp=lp,&
1132 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1133 0 : ELSE IF(ALLOCATED(this%uu)) THEN
1134 : this%uu(iz,:,:,:,ipm) = rotMMPmat(this%uu(iz,:,:,:,ipm),alpha,beta,gamma,l,lp=lp,&
1135 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1136 : this%ud(iz,:,:,:,ipm) = rotMMPmat(this%ud(iz,:,:,:,ipm),alpha,beta,gamma,l,lp=lp,&
1137 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1138 : this%du(iz,:,:,:,ipm) = rotMMPmat(this%du(iz,:,:,:,ipm),alpha,beta,gamma,l,lp=lp,&
1139 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1140 : this%dd(iz,:,:,:,ipm) = rotMMPmat(this%dd(iz,:,:,:,ipm),alpha,beta,gamma,l,lp=lp,&
1141 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1142 :
1143 0 : IF(ALLOCATED(this%uulo)) THEN
1144 0 : iLO_ind = 0
1145 0 : DO ilo = 1, atoms%nlo(atomType)
1146 0 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
1147 0 : iLO_ind = iLO_ind + 1
1148 : this%uulo(iz,:,:,iLO_ind,:,ipm) = rotMMPmat(this%uulo(iz,:,:,iLO_ind,:,ipm),alpha,beta,gamma,l,lp=lp,&
1149 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1150 : this%dulo(iz,:,:,iLO_ind,:,ipm) = rotMMPmat(this%dulo(iz,:,:,iLO_ind,:,ipm),alpha,beta,gamma,l,lp=lp,&
1151 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1152 : ENDDO
1153 0 : iLO_ind = 0
1154 0 : DO ilo = 1, atoms%nlo(atomTypep)
1155 0 : IF(atoms%llo(ilo,atomTypep).NE.lp) CYCLE
1156 0 : iLO_ind = iLO_ind + 1
1157 : this%ulou(iz,:,:,iLO_ind,:,ipm) = rotMMPmat(this%ulou(iz,:,:,iLO_ind,:,ipm),alpha,beta,gamma,l,lp=lp,&
1158 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1159 : this%ulod(iz,:,:,iLO_ind,:,ipm) = rotMMPmat(this%ulod(iz,:,:,iLO_ind,:,ipm),alpha,beta,gamma,l,lp=lp,&
1160 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1161 : ENDDO
1162 0 : iLO_ind = 0
1163 0 : DO ilo = 1, atoms%nlo(atomType)
1164 0 : IF(atoms%llo(ilo,atomType).NE.l) CYCLE
1165 0 : iLOp_ind = 0
1166 0 : iLO_ind = iLO_ind + 1
1167 0 : DO ilop = 1, atoms%nlo(atomTypep)
1168 0 : IF(atoms%llo(ilop,atomTypep).NE.lp) CYCLE
1169 0 : iLOp_ind = iLOp_ind + 1
1170 : this%uloulop(iz,:,:,iLO_ind,iLOp_ind,:,ipm) = rotMMPmat(this%uloulop(iz,:,:,iLO_ind,iLOp_ind,:,ipm),alpha,beta,gamma,l,lp=lp,&
1171 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1172 : ENDDO
1173 : ENDDO
1174 : ENDIF
1175 0 : ELSE IF(ALLOCATED(this%gmmpMat_k)) THEN
1176 0 : DO ikpt = 1, SIZE(this%gmmpMat_k,6)
1177 : this%gmmpMat_k(iz,:,:,:,ipm,ikpt) = rotMMPmat(this%gmmpMat_k(iz,:,:,:,ipm,ikpt),alpha,beta,gamma,l,lp=lp,&
1178 0 : spin_rotation=spin_rotation,real_space_rotation=real_space_rotation)
1179 : ENDDO
1180 : ENDIF
1181 : ENDDO
1182 : ENDDO
1183 0 : CALL timestop("Green's Function: Rotate (Angles)")
1184 : END SUBROUTINE rotate_euler_angles_gf
1185 :
1186 0 : PURE FUNCTION checkEmpty_greensf(this,m,mp,spin,ipm) Result(l_empty)
1187 :
1188 : CLASS(t_greensf), INTENT(IN) :: this
1189 : INTEGER, INTENT(IN) :: m
1190 : INTEGER, INTENT(IN) :: mp
1191 : INTEGER, INTENT(IN) :: spin
1192 : INTEGER, INTENT(IN) :: ipm
1193 :
1194 : LOGICAL :: l_empty
1195 :
1196 0 : IF(ALLOCATED(this%gmmpMat)) THEN
1197 0 : l_empty = ALL(ABS(this%gmmpMat(:,m,mp,spin,ipm)).LT.1e-12)
1198 0 : ELSE IF(ALLOCATED(this%uu)) THEN
1199 : l_empty = ALL(ABS(this%uu(:,m,mp,spin,ipm)).LT.1e-12) &
1200 : .AND. ALL(ABS(this%dd(:,m,mp,spin,ipm)).LT.1e-12) &
1201 : .AND. ALL(ABS(this%ud(:,m,mp,spin,ipm)).LT.1e-12) &
1202 0 : .AND. ALL(ABS(this%du(:,m,mp,spin,ipm)).LT.1e-12)
1203 0 : IF(ALLOCATED(this%uulo)) THEN
1204 : l_empty = l_empty .AND. ALL(ABS(this%uulo(:,m,mp,:,spin,ipm)).LT.1e-12) &
1205 : .AND. ALL(ABS(this%ulou(:,m,mp,:,spin,ipm)).LT.1e-12) &
1206 : .AND. ALL(ABS(this%ulod(:,m,mp,:,spin,ipm)).LT.1e-12) &
1207 : .AND. ALL(ABS(this%dulo(:,m,mp,:,spin,ipm)).LT.1e-12) &
1208 0 : .AND. ALL(ABS(this%uloulop(:,m,mp,:,:,spin,ipm)).LT.1e-12)
1209 : ENDIF
1210 0 : ELSE IF(ALLOCATED(this%gmmpMat_k)) THEN
1211 0 : l_empty = ALL(ABS(this%gmmpMat_k(:,m,mp,spin,ipm,:)).LT.1e-12)
1212 : ENDIF
1213 :
1214 0 : END FUNCTION checkEmpty_greensf
1215 :
1216 1018 : SUBROUTINE reset_gf(this)
1217 :
1218 : !---------------------------------------------------
1219 : ! Sets all gmmpMat arrays back to 0
1220 : !---------------------------------------------------
1221 :
1222 : CLASS(t_greensf), INTENT(INOUT) :: this
1223 :
1224 27656454 : IF(ALLOCATED(this%gmmpMat)) this%gmmpMat = cmplx_0
1225 1018 : IF(ALLOCATED(this%uu)) THEN
1226 303828 : this%uu = cmplx_0
1227 303828 : this%ud = cmplx_0
1228 303828 : this%du = cmplx_0
1229 303828 : this%dd = cmplx_0
1230 : ENDIF
1231 1018 : IF(ALLOCATED(this%uulo)) THEN
1232 253216 : this%uulo = cmplx_0
1233 253216 : this%ulou = cmplx_0
1234 253216 : this%dulo = cmplx_0
1235 253216 : this%ulod = cmplx_0
1236 :
1237 354520 : this%uloulop = cmplx_0
1238 : ENDIF
1239 1018 : IF(ALLOCATED(this%gmmpMat_k)) this%gmmpMat_k = cmplx_0
1240 :
1241 1018 : END SUBROUTINE reset_gf
1242 :
1243 :
1244 0 : SUBROUTINE resetSingleElem_gf(this,m,mp,spin,ipm)
1245 :
1246 : !---------------------------------------------------
1247 : ! Sets one Element in gmmpMat arrays back to 0
1248 : !---------------------------------------------------
1249 :
1250 : CLASS(t_greensf), INTENT(INOUT) :: this
1251 : INTEGER, INTENT(IN) :: m
1252 : INTEGER, INTENT(IN) :: mp
1253 : INTEGER, INTENT(IN) :: spin
1254 : INTEGER, INTENT(IN) :: ipm
1255 :
1256 0 : IF(ALLOCATED(this%gmmpMat)) this%gmmpMat(:,m,mp,spin,ipm) = cmplx_0
1257 0 : IF(ALLOCATED(this%uu)) THEN
1258 0 : this%uu(:,m,mp,spin,ipm) = cmplx_0
1259 0 : this%ud(:,m,mp,spin,ipm) = cmplx_0
1260 0 : this%du(:,m,mp,spin,ipm) = cmplx_0
1261 0 : this%dd(:,m,mp,spin,ipm) = cmplx_0
1262 : ENDIF
1263 0 : IF(ALLOCATED(this%uulo)) THEN
1264 0 : this%uulo(:,m,mp,:,spin,ipm) = cmplx_0
1265 0 : this%ulou(:,m,mp,:,spin,ipm) = cmplx_0
1266 0 : this%dulo(:,m,mp,:,spin,ipm) = cmplx_0
1267 0 : this%ulod(:,m,mp,:,spin,ipm) = cmplx_0
1268 :
1269 0 : this%uloulop(:,m,mp,:,:,spin,ipm) = cmplx_0
1270 : ENDIF
1271 0 : IF(ALLOCATED(this%gmmpMat_k)) this%gmmpMat_k(:,m,mp,spin,ipm,:) = cmplx_0
1272 :
1273 0 : END SUBROUTINE resetSingleElem_gf
1274 :
1275 0 : FUNCTION integrateOverMT_greensf(this,atoms,input,gfinp,f,g,flo,l_fullRadial) Result(gIntegrated)
1276 :
1277 : USE m_intgr
1278 : USE m_types_usdus
1279 : USE m_types_denCoeffsOffDiag
1280 : USE m_types_mat
1281 :
1282 : CLASS(t_greensf), INTENT(IN) :: this
1283 : TYPE(t_atoms), INTENT(IN) :: atoms
1284 : TYPE(t_input), INTENT(IN) :: input
1285 : TYPE(t_gfinp), INTENT(IN) :: gfinp
1286 : REAL, INTENT(IN) :: f(:,:,0:,:,:)
1287 : REAL, INTENT(IN) :: g(:,:,0:,:,:)
1288 : REAL, INTENT(IN) :: flo(:,:,:,:,:)
1289 : LOGICAL, OPTIONAL, INTENT(IN) :: l_fullRadial
1290 :
1291 : TYPE(t_greensf) :: gIntegrated
1292 :
1293 : LOGICAL :: l_fullRadialArg,l_explicit
1294 : INTEGER :: l,lp,atomType,atomTypep,ipm,spin,m,mp,iz,jr,jrp
1295 : REAL :: realPart, imagPart
1296 0 : COMPLEX, ALLOCATABLE :: gmatR(:,:)
1297 0 : COMPLEX :: gmat(atoms%jmtd)
1298 0 : TYPE(t_mat) :: gmatTmp
1299 :
1300 0 : l_fullRadialArg = .FALSE.
1301 0 : IF(PRESENT(l_fullRadial)) l_fullRadialArg = l_fullRadial
1302 :
1303 0 : IF(this%elem%l_sphavg) CALL juDFT_error("GF has to be provided with radial dependence",&
1304 0 : calledby="integrateOverMT_greensf")
1305 :
1306 :
1307 0 : CALL timestart("Green's Function: Average over MT")
1308 0 : CALL gIntegrated%init(this%elem,gfinp,atoms,input,contour_in=this%contour,l_sphavg_in=.TRUE.)
1309 0 : gIntegrated%l_calc = .TRUE.
1310 0 : l = this%elem%l
1311 0 : lp = this%elem%lp
1312 0 : atomType = this%elem%atomType
1313 0 : atomTypep = this%elem%atomTypep
1314 : !Do we have the offdiagonal scalar products
1315 0 : l_explicit = .TRUE.
1316 0 : IF(ALLOCATED(this%scalarProducts%uun)) l_explicit = .FALSE.
1317 :
1318 : !only intersite arguments have independent radial arguments ??
1319 0 : l_fullRadialArg = l_fullRadialArg.AND.this%elem%isIntersite()
1320 :
1321 0 : DO ipm = 1, 2
1322 0 : DO spin = 1 , SIZE(this%uu,4)
1323 0 : IF(.NOT.l_explicit) THEN
1324 0 : DO iz = 1, this%contour%nz
1325 0 : CALL this%get(atoms,iz,ipm==2,spin,gmatTmp)
1326 0 : CALL gIntegrated%set(iz,ipm==2,gmatTmp,spin=spin)
1327 : ENDDO
1328 : ELSE
1329 0 : DO mp = -lp, lp
1330 0 : DO m = -l, l
1331 0 : IF(this%checkEmpty(m,mp,spin,ipm)) CYCLE
1332 0 : IF(.NOT.l_fullRadialArg) THEN
1333 0 : CALL this%getRadial(atoms,m,mp,ipm==2,spin,f,g,flo,gmatR)
1334 : ENDIF
1335 0 : DO iz = 1, this%contour%nz
1336 0 : IF(l_fullRadialArg) THEN
1337 0 : CALL this%getRadialRadial(atoms,iz,m,mp,ipm==2,spin,f,g,flo,gmatR)
1338 0 : gmat = cmplx_0
1339 0 : DO jr = 1, SIZE(gmat)
1340 0 : CALL intgr3(REAL(gmatR(:,jr)),atoms%rmsh(:,atomTypep),atoms%dx(atomTypep),atoms%jri(atomTypep),realPart)
1341 0 : CALL intgr3(AIMAG(gmatR(:,jr)),atoms%rmsh(:,atomTypep),atoms%dx(atomTypep),atoms%jri(atomTypep),imagPart)
1342 :
1343 0 : gmat(jr) = realPart + ImagUnit * imagPart
1344 : ENDDO
1345 : ELSE
1346 0 : gmat = gmatR(:,iz)
1347 : ENDIF
1348 0 : CALL intgr3(REAL(gmat),atoms%rmsh(:,atomType),atoms%dx(atomType),atoms%jri(atomType),realPart)
1349 0 : CALL intgr3(AIMAG(gmat),atoms%rmsh(:,atomType),atoms%dx(atomType),atoms%jri(atomType),imagPart)
1350 :
1351 0 : gIntegrated%gmmpMat(iz,m,mp,spin,ipm) = realPart + ImagUnit * imagPart
1352 : ENDDO
1353 : ENDDO
1354 : ENDDO
1355 : ENDIF
1356 : ENDDO
1357 : ENDDO
1358 0 : CALL timestop("Green's Function: Average over MT")
1359 :
1360 0 : END FUNCTION integrateOverMT_greensf
1361 :
1362 0 : END MODULE m_types_greensf
|
