Line data Source code
1 : MODULE m_types_greensfContourData
2 : USE m_juDFT
3 : USE m_types_gfinp
4 : USE m_constants
5 : IMPLICIT NONE
6 : PRIVATE
7 :
8 : TYPE t_greensfContourData
9 : !This type is used in types_greensf for the storage of the energy points and weights
10 : !It is defined here, because we want to use it in eContour_gfinp, without importing types_greensf
11 : INTEGER :: nz = 0 !Number of points
12 : COMPLEX, ALLOCATABLE :: e(:) !energy points
13 : COMPLEX, ALLOCATABLE :: de(:) !integration weights
14 :
15 : CONTAINS
16 : PROCEDURE,PASS :: init => init_greensfContourData
17 : PROCEDURE :: eContour => eContour_greensfContourData
18 : PROCEDURE :: mpi_bc => mpi_bc_greensfContourData
19 : END TYPE t_greensfContourData
20 :
21 : PUBLIC t_greensfContourData
22 :
23 : CONTAINS
24 :
25 1058 : SUBROUTINE init_greensfContourData(this,contourInp,contour_in)
26 :
27 : CLASS(t_greensfContourData), INTENT(INOUT) :: this
28 : TYPE(t_contourInp), INTENT(IN) :: contourInp
29 : TYPE(t_greensfContourData), OPTIONAL, INTENT(IN) :: contour_in
30 : !
31 : !Setting up parameters for the energy contour if one was passed
32 : !
33 1058 : IF(PRESENT(contour_in)) THEN
34 0 : this%nz = contour_in%nz
35 0 : ALLOCATE(this%e(this%nz),source=cmplx_0)
36 0 : ALLOCATE(this%de(this%nz),source=cmplx_0)
37 :
38 0 : this%e = contour_in%e
39 0 : this%de = contour_in%de
40 : ELSE
41 1062 : SELECT CASE(contourInp%shape)
42 : CASE(CONTOUR_RECTANGLE_CONST)
43 4 : this%nz = contourInp%n1+contourInp%n2+contourInp%n3+contourInp%nmatsub
44 : CASE(CONTOUR_SEMICIRCLE_CONST)
45 1050 : this%nz = contourInp%ncirc
46 : CASE(CONTOUR_DOS_CONST)
47 4 : this%nz = contourInp%nDOS
48 : CASE DEFAULT
49 1058 : CALL juDFT_error("No valid energy contour mode",calledby="init_contourDataType")
50 : END SELECT
51 :
52 159866 : ALLOCATE(this%e(this%nz),source=cmplx_0)
53 158808 : ALLOCATE(this%de(this%nz),source=cmplx_0)
54 : ENDIF
55 :
56 1058 : END SUBROUTINE init_greensfContourData
57 :
58 1018 : SUBROUTINE mpi_bc_greensfContourData(this,mpi_comm,irank)
59 : USE m_mpi_bc_tool
60 : CLASS(t_greensfContourData), INTENT(INOUT)::this
61 : INTEGER, INTENT(IN):: mpi_comm
62 : INTEGER, INTENT(IN), OPTIONAL::irank
63 : INTEGER ::rank
64 1018 : IF (PRESENT(irank)) THEN
65 1018 : rank = irank
66 : ELSE
67 0 : rank = 0
68 : END IF
69 :
70 1018 : CALL mpi_bc(this%nz,rank,mpi_comm)
71 1018 : CALL mpi_bc(this%e,rank,mpi_comm)
72 1018 : CALL mpi_bc(this%de,rank,mpi_comm)
73 :
74 1018 : END SUBROUTINE mpi_bc_greensfContourData
75 :
76 46 : SUBROUTINE eContour_greensfContourData(this,contourInp,ef,irank)
77 :
78 : USE m_grule
79 :
80 : !Calculates the complex energy contour and
81 : !writes it into the corresponding arrays in gf
82 :
83 : CLASS(t_greensfContourData), INTENT(INOUT) :: this
84 : TYPE(t_contourInp), INTENT(IN) :: contourInp
85 : REAL, INTENT(IN) :: ef
86 : INTEGER, INTENT(IN) :: irank
87 :
88 : INTEGER iz,n
89 : REAL e1, e2, sigma
90 : COMPLEX del
91 : REAL r, xr, expo, ff
92 46 : REAL, ALLOCATABLE :: x(:), w(:)
93 :
94 : !Help arrays
95 16660 : ALLOCATE(x(this%nz),source=0.0)
96 16614 : ALLOCATE(w(this%nz),source=0.0)
97 :
98 :
99 : !Transform from relative to ef to absolute
100 46 : e1 = ef+contourInp%eb
101 46 : e2 = ef+contourInp%et
102 :
103 48 : SELECT CASE(contourInp%shape)
104 :
105 : CASE(CONTOUR_RECTANGLE_CONST)
106 2 : sigma = contourInp%sigma * pi_const
107 2 : IF(contourInp%nmatsub > 0) THEN
108 2 : n = 0
109 :
110 : !Left Vertical part (e1,0) -> (e1,sigma)
111 2 : del = contourInp%nmatsub * CMPLX(0.0,sigma)
112 2 : CALL grule(contourInp%n1,x(1:(contourInp%n1)/2),w(1:(contourInp%n1)/2))
113 348 : x = -x
114 12 : DO iz = 1, (contourInp%n1+3)/2-1
115 10 : x(contourInp%n1-iz+1) = -x(iz)
116 12 : w(contourInp%n1-iz+1) = w(iz)
117 : ENDDO
118 22 : DO iz = 1, contourInp%n1
119 20 : n = n + 1
120 20 : IF(n.GT.this%nz) CALL juDFT_error("Dimension error in energy mesh",calledby="eContour_gfinp")
121 20 : this%e(n) = e1 + del + del * x(iz)
122 22 : this%de(n) = w(iz)*del
123 : ENDDO
124 :
125 : !Horizontal Part (eb,sigma) -> (et,sigma)
126 2 : del = (ef-30*contourInp%sigma-e1)/2.0
127 2 : CALL grule(contourInp%n2,x(1:(contourInp%n2)/2),w(1:(contourInp%n2)/2))
128 348 : x = -x
129 130 : DO iz = 1, (contourInp%n2+3)/2-1
130 128 : x(contourInp%n2-iz+1) = -x(iz)
131 130 : w(contourInp%n2-iz+1) = w(iz)
132 : ENDDO
133 258 : DO iz = 1, contourInp%n2
134 256 : n = n + 1
135 256 : IF(n.GT.this%nz) CALL juDFT_error("Dimension error in energy mesh",calledby="eContour_gfinp")
136 256 : this%e(n) = del*x(iz) + del + e1 + 2 * contourInp%nmatsub * ImagUnit * sigma
137 258 : this%de(n) = del*w(iz)
138 : ENDDO
139 :
140 : !Right Vertical part (et,sigma) -> infty
141 2 : CALL grule(contourInp%n3,x(1:(contourInp%n3)/2),w(1:(contourInp%n3)/2))
142 348 : x = -x
143 32 : DO iz = 1, (contourInp%n3+3)/2-1
144 30 : x(contourInp%n3-iz+1) = -x(iz)
145 32 : w(contourInp%n3-iz+1) = w(iz)
146 : ENDDO
147 2 : del = 30*contourInp%sigma
148 62 : DO iz = 1, contourInp%n3
149 60 : n = n + 1
150 60 : IF(n.GT.this%nz) CALL juDFT_error("Dimension error in energy mesh",calledby="eContour_gfinp")
151 60 : this%e(n) = del*x(iz)+ef + 2 * contourInp%nmatsub * ImagUnit * sigma
152 60 : expo = -ABS(REAL(this%e(n))-ef)/contourInp%sigma
153 60 : expo = EXP(expo)
154 60 : IF(REAL(this%e(n))<ef) THEN
155 30 : ff = 1.0/(expo+1.0)
156 : ELSE
157 30 : ff = expo/(expo+1.0)
158 : ENDIF
159 62 : this%de(n) = w(iz)*del * ff
160 : ENDDO
161 :
162 : !Matsubara frequencies
163 12 : DO iz = contourInp%nmatsub , 1, -1
164 10 : n = n + 1
165 10 : IF(n.GT.this%nz) CALL juDFT_error("Dimension error in energy mesh",calledby="eContour_gfinp")
166 10 : this%e(n) = ef + (2*iz-1) * ImagUnit *sigma
167 12 : this%de(n) = -2 * ImagUnit * sigma
168 : ENDDO
169 : ENDIF
170 : CASE(CONTOUR_SEMICIRCLE_CONST)
171 :
172 : !Radius
173 42 : r = (e2-e1)*0.5
174 : !midpoint
175 42 : xr = (e2+e1)*0.5
176 :
177 42 : CALL grule(contourInp%ncirc,x(1:(contourInp%ncirc)/2),w(1:(contourInp%ncirc)/2))
178 :
179 2730 : DO iz = 1, contourInp%ncirc/2
180 2688 : x(contourInp%ncirc-iz+1) = -x(iz)
181 2730 : w(contourInp%ncirc-iz+1) = w(iz)
182 : ENDDO
183 5418 : DO iz = 1, contourInp%ncirc
184 5376 : this%e(iz) = xr + ImagUnit * r * EXP(ImagUnit*pi_const/2.0 * x(iz))
185 5376 : this%de(iz) = pi_const/2.0 * r * w(iz) * EXP(ImagUnit*pi_const/2.0 * x(iz))
186 : !Scale the imaginary part with the given factor alpha
187 5376 : this%e(iz) = REAL(this%e(iz)) + ImagUnit * contourInp%alpha * AIMAG(this%e(iz))
188 5418 : this%de(iz) = REAL(this%de(iz)) + ImagUnit * contourInp%alpha * AIMAG(this%de(iz))
189 : ENDDO
190 :
191 : CASE(CONTOUR_DOS_CONST)
192 :
193 : !Equidistant contour (without vertical edges)
194 2 : del = (contourInp%et-contourInp%eb)/REAL(this%nz-1)
195 10802 : DO iz = 1, this%nz
196 10800 : this%e(iz) = (iz-1) * del + e1 + ImagUnit * contourInp%sigmaDOS
197 10802 : IF(contourInp%l_dosfermi) THEN
198 0 : expo = -ABS(REAL(this%e(iz))-ef)/contourInp%sigmaDOS
199 0 : expo = EXP(expo)
200 0 : IF(REAL(this%e(iz))<ef) THEN
201 0 : ff = 1.0/(expo+1.0)
202 : ELSE
203 0 : ff = expo/(expo+1.0)
204 : ENDIF
205 0 : this%de(iz) = del * ff
206 : ELSE
207 10800 : this%de(iz) = del
208 : ENDIF
209 : ENDDO
210 :
211 : !Not really important but for trapezian method
212 : !the weight is half at the edges
213 2 : this%de(1) = this%de(1)/2.0
214 2 : this%de(this%nz) = this%de(this%nz)/2.0
215 :
216 : CASE DEFAULT
217 46 : CALL juDFT_error("Invalid mode for energy contour in Green's function calculation", calledby="eContour_gfinp")
218 : END SELECT
219 :
220 46 : IF(irank.EQ.0) THEN
221 : !Write out the information about the energy contour
222 23 : WRITE(oUnit,"(A)") "---------------------------------------------"
223 23 : WRITE(oUnit,"(A)") " Green's function energy contour"
224 23 : WRITE(oUnit,"(A)") "---------------------------------------------"
225 23 : WRITE(oUnit,999) TRIM(ADJUSTL(contourInp%label))
226 23 : WRITE(oUnit,1000) contourInp%shape
227 :
228 1 : SELECT CASE(contourInp%shape)
229 :
230 : CASE(CONTOUR_RECTANGLE_CONST)
231 1 : WRITE(oUnit,"(A)") "Rectangular Contour: "
232 1 : WRITE(oUnit,1010) this%nz, contourInp%nmatsub,contourInp%n1,contourInp%n2,contourInp%n3
233 1 : WRITE(oUnit,"(A)") "Energy limits (rel. to fermi energy): "
234 1 : WRITE(oUnit,1040) contourInp%eb,0.0
235 : CASE(CONTOUR_SEMICIRCLE_CONST)
236 21 : WRITE(oUnit,"(A)") "Semicircle Contour: "
237 21 : WRITE(oUnit,1020) this%nz, contourInp%alpha
238 21 : WRITE(oUnit,"(A)") "Energy limits (rel. to fermi energy): "
239 21 : WRITE(oUnit,1040) contourInp%eb,contourInp%et
240 : CASE(CONTOUR_DOS_CONST)
241 1 : WRITE(oUnit,"(A)") "Equidistant Contour for DOS calculations: "
242 1 : WRITE(oUnit,1030) this%nz, contourInp%sigmaDOS
243 1 : WRITE(oUnit,"(A)") "Energy limits (rel. to fermi energy): "
244 24 : WRITE(oUnit,1040) contourInp%eb,contourInp%et
245 : CASE DEFAULT
246 :
247 : END SELECT
248 :
249 : !Write out points and weights
250 23 : WRITE(oUnit,*)
251 23 : WRITE(oUnit,"(A)") " Energy points: "
252 23 : WRITE(oUnit,"(A)") "---------------------------------------------"
253 8284 : DO iz = 1, this%nz
254 8284 : WRITE(oUnit,1050) REAL(this%e(iz)), AIMAG(this%e(iz)), REAL(this%de(iz)), AIMAG(this%de(iz))
255 : ENDDO
256 :
257 : 999 FORMAT("Name of the contour:", A)
258 : 1000 FORMAT("Using energy contour mode: ", I1,/)
259 : 1010 FORMAT("nz: ", I5.1,"; nmatsub: ", I5.1,"; n1: ", I5.1,"; n2: ", I5.1,"; n3: ", I5.1)
260 : 1020 FORMAT("nz: ", I5.1," alpha: ", f8.4)
261 : 1030 FORMAT("n: ", I5.1,"; sigma: ", f8.4)
262 : 1040 FORMAT("eb: ", f8.4,"; et: ",f8.4)
263 : 1050 FORMAT(2f8.4," weight: ",2e15.4)
264 : ENDIF
265 :
266 46 : END SUBROUTINE eContour_greensfContourData
267 :
268 :
269 0 : END MODULE m_types_greensfContourData
|
