Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Subroutines to compute Green functions
10 : ! **************************************************************************************************
11 : MODULE negf_green_methods
12 : USE cp_cfm_basic_linalg, ONLY: cp_cfm_gemm,&
13 : cp_cfm_lu_invert,&
14 : cp_cfm_norm,&
15 : cp_cfm_scale,&
16 : cp_cfm_scale_and_add,&
17 : cp_cfm_scale_and_add_fm,&
18 : cp_cfm_transpose
19 : USE cp_cfm_types, ONLY: cp_cfm_create,&
20 : cp_cfm_get_info,&
21 : cp_cfm_release,&
22 : cp_cfm_to_cfm,&
23 : cp_cfm_type,&
24 : cp_fm_to_cfm
25 : USE cp_fm_struct, ONLY: cp_fm_struct_get,&
26 : cp_fm_struct_type
27 : USE cp_fm_types, ONLY: cp_fm_get_info,&
28 : cp_fm_type
29 : USE kinds, ONLY: dp
30 : USE mathconstants, ONLY: gaussi,&
31 : z_mone,&
32 : z_one,&
33 : z_zero
34 : USE parallel_gemm_api, ONLY: parallel_gemm
35 : #include "./base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 : PRIVATE
39 :
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'negf_green_methods'
41 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
42 :
43 : PUBLIC :: sancho_work_matrices_type, sancho_work_matrices_create, sancho_work_matrices_release
44 : PUBLIC :: do_sancho, negf_contact_self_energy, negf_contact_broadening_matrix, negf_retarded_green_function
45 :
46 : TYPE sancho_work_matrices_type
47 : ! A_{n+1} = A_n + D_n + E_n
48 : TYPE(cp_cfm_type), POINTER :: a => NULL()
49 : ! A0_{n+1} = A0_n + D_n
50 : TYPE(cp_cfm_type), POINTER :: a0 => NULL()
51 : ! A_inv = A^-1
52 : TYPE(cp_cfm_type), POINTER :: a_inv => NULL()
53 : ! B_{n+1} = - B_n A_n^-1 B_n \equiv B_n A_n^-1 B_n
54 : TYPE(cp_cfm_type), POINTER :: b => NULL()
55 : ! C_{n+1} = - C_n A_n^-1 C_n \equiv C_n A_n^-1 C_n
56 : TYPE(cp_cfm_type), POINTER :: c => NULL()
57 : ! D_n = - B_n A_n^-1 C_n
58 : TYPE(cp_cfm_type), POINTER :: d => NULL()
59 : ! E_n = - C_n A_n^-1 B_n
60 : TYPE(cp_cfm_type), POINTER :: e => NULL()
61 : ! a scratch area for matrix multiplication
62 : TYPE(cp_cfm_type), POINTER :: scratch => NULL()
63 : END TYPE sancho_work_matrices_type
64 :
65 : CONTAINS
66 : ! **************************************************************************************************
67 : !> \brief Create work matrices required for the Lopez-Sancho algorithm.
68 : !> \param work work matrices to create (allocated and initialised on exit)
69 : !> \param fm_struct dense matrix structure
70 : !> \author Sergey Chulkov
71 : ! **************************************************************************************************
72 2448 : SUBROUTINE sancho_work_matrices_create(work, fm_struct)
73 : TYPE(sancho_work_matrices_type), INTENT(inout) :: work
74 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
75 :
76 : CHARACTER(len=*), PARAMETER :: routineN = 'sancho_work_matrices_create'
77 :
78 : INTEGER :: handle, ncols, nrows
79 :
80 1224 : CALL timeset(routineN, handle)
81 1224 : CPASSERT(ASSOCIATED(fm_struct))
82 :
83 1224 : CALL cp_fm_struct_get(fm_struct, nrow_global=nrows, ncol_global=ncols)
84 1224 : CPASSERT(nrows == ncols)
85 :
86 1224 : NULLIFY (work%a, work%a0, work%a_inv, work%b, work%c, work%d, work%e, work%scratch)
87 1224 : ALLOCATE (work%a, work%a0, work%a_inv, work%b, work%c, work%d, work%e, work%scratch)
88 1224 : CALL cp_cfm_create(work%a, fm_struct)
89 1224 : CALL cp_cfm_create(work%a0, fm_struct)
90 1224 : CALL cp_cfm_create(work%a_inv, fm_struct)
91 1224 : CALL cp_cfm_create(work%b, fm_struct)
92 1224 : CALL cp_cfm_create(work%c, fm_struct)
93 1224 : CALL cp_cfm_create(work%d, fm_struct)
94 1224 : CALL cp_cfm_create(work%e, fm_struct)
95 1224 : CALL cp_cfm_create(work%scratch, fm_struct)
96 :
97 1224 : CALL timestop(handle)
98 1224 : END SUBROUTINE sancho_work_matrices_create
99 :
100 : ! **************************************************************************************************
101 : !> \brief Release work matrices.
102 : !> \param work work matrices to release
103 : !> \author Sergey Chulkov
104 : ! **************************************************************************************************
105 1224 : SUBROUTINE sancho_work_matrices_release(work)
106 : TYPE(sancho_work_matrices_type), INTENT(inout) :: work
107 :
108 : CHARACTER(len=*), PARAMETER :: routineN = 'sancho_work_matrices_release'
109 :
110 : INTEGER :: handle
111 :
112 1224 : CALL timeset(routineN, handle)
113 :
114 1224 : IF (ASSOCIATED(work%a)) THEN
115 1224 : CALL cp_cfm_release(work%a)
116 1224 : DEALLOCATE (work%a)
117 : NULLIFY (work%a)
118 : END IF
119 1224 : IF (ASSOCIATED(work%a0)) THEN
120 1224 : CALL cp_cfm_release(work%a0)
121 1224 : DEALLOCATE (work%a0)
122 : NULLIFY (work%a0)
123 : END IF
124 1224 : IF (ASSOCIATED(work%a_inv)) THEN
125 1224 : CALL cp_cfm_release(work%a_inv)
126 1224 : DEALLOCATE (work%a_inv)
127 : NULLIFY (work%a_inv)
128 : END IF
129 1224 : IF (ASSOCIATED(work%b)) THEN
130 1224 : CALL cp_cfm_release(work%b)
131 1224 : DEALLOCATE (work%b)
132 : NULLIFY (work%b)
133 : END IF
134 1224 : IF (ASSOCIATED(work%c)) THEN
135 1224 : CALL cp_cfm_release(work%c)
136 1224 : DEALLOCATE (work%c)
137 : NULLIFY (work%c)
138 : END IF
139 1224 : IF (ASSOCIATED(work%d)) THEN
140 1224 : CALL cp_cfm_release(work%d)
141 1224 : DEALLOCATE (work%d)
142 : NULLIFY (work%d)
143 : END IF
144 1224 : IF (ASSOCIATED(work%e)) THEN
145 1224 : CALL cp_cfm_release(work%e)
146 1224 : DEALLOCATE (work%e)
147 : NULLIFY (work%e)
148 : END IF
149 1224 : IF (ASSOCIATED(work%scratch)) THEN
150 1224 : CALL cp_cfm_release(work%scratch)
151 1224 : DEALLOCATE (work%scratch)
152 : NULLIFY (work%scratch)
153 : END IF
154 :
155 1224 : CALL timestop(handle)
156 1224 : END SUBROUTINE sancho_work_matrices_release
157 :
158 : ! **************************************************************************************************
159 : !> \brief Iterative method to compute a retarded surface Green's function at the point omega.
160 : !> \param g_surf computed retarded surface Green's function (initialised on exit)
161 : !> \param omega argument of the Green's function
162 : !> \param h0 diagonal block of the Kohn-Sham matrix (must be Hermitian)
163 : !> \param s0 diagonal block of the overlap matrix (must be Hermitian)
164 : !> \param h1 off-fiagonal block of the Kohn-Sham matrix
165 : !> \param s1 off-fiagonal block of the overlap matrix
166 : !> \param conv convergence threshold
167 : !> \param transp flag which indicates that the matrices h1 and s1 matrices should be transposed
168 : !> \param work a set of work matrices needed to compute the surface Green's function
169 : !> \author Sergey Chulkov
170 : ! **************************************************************************************************
171 13220 : SUBROUTINE do_sancho(g_surf, omega, h0, s0, h1, s1, conv, transp, work)
172 : TYPE(cp_cfm_type), INTENT(IN) :: g_surf
173 : COMPLEX(kind=dp), INTENT(in) :: omega
174 : TYPE(cp_fm_type), INTENT(IN) :: h0, s0, h1, s1
175 : REAL(kind=dp), INTENT(in) :: conv
176 : LOGICAL, INTENT(in) :: transp
177 : TYPE(sancho_work_matrices_type), INTENT(in) :: work
178 :
179 : CHARACTER(len=*), PARAMETER :: routineN = 'do_sancho'
180 :
181 : INTEGER :: handle, ncols, nrows
182 :
183 6610 : CALL timeset(routineN, handle)
184 :
185 6610 : CALL cp_cfm_get_info(g_surf, nrow_global=nrows, ncol_global=ncols)
186 :
187 : IF (debug_this_module) THEN
188 6610 : CPASSERT(nrows == ncols)
189 : END IF
190 :
191 : ! A^{ret.}_0 = omega * s0 - h0
192 6610 : CALL cp_fm_to_cfm(msourcer=s0, mtarget=work%a)
193 6610 : CALL cp_cfm_scale_and_add_fm(omega, work%a, z_mone, h0)
194 : ! A0^{ret.}_0 = A^{ret.}_0
195 6610 : CALL cp_cfm_to_cfm(work%a, work%a0)
196 :
197 : ! B^{ret.}_0 = omega * s1 - h1
198 : ! C^{ret.}_0 = B^{ret.}_0^T
199 6610 : IF (transp) THEN
200 : ! beta = omega * s1 - h1
201 1260 : CALL cp_fm_to_cfm(msourcer=s1, mtarget=work%c)
202 1260 : CALL cp_cfm_scale_and_add_fm(omega, work%c, z_mone, h1)
203 :
204 : ! alpha = omega * s1' - h1'
205 1260 : CALL cp_cfm_transpose(matrix=work%c, trans='T', matrixt=work%b)
206 : ELSE
207 : ! alpha = omega * s1 - h1
208 5350 : CALL cp_fm_to_cfm(msourcer=s1, mtarget=work%b)
209 5350 : CALL cp_cfm_scale_and_add_fm(omega, work%b, z_mone, h1)
210 :
211 : ! beta = omega * s1' - h1'
212 5350 : CALL cp_cfm_transpose(matrix=work%b, trans='T', matrixt=work%c)
213 : END IF
214 :
215 : ! actually compute the Green's function
216 44026 : DO WHILE (cp_cfm_norm(work%b, 'M') + cp_cfm_norm(work%c, 'M') > conv)
217 : ! A_n^-1
218 37416 : CALL cp_cfm_to_cfm(work%a, work%a_inv)
219 37416 : CALL cp_cfm_lu_invert(work%a_inv)
220 :
221 : ! scratch <- A_n^-1 * B_n
222 37416 : CALL parallel_gemm('N', 'N', nrows, nrows, nrows, z_one, work%a_inv, work%b, z_zero, work%scratch)
223 : ! E_n = - C_n A_n^-1 B_n
224 37416 : CALL cp_cfm_gemm('N', 'N', nrows, nrows, nrows, z_mone, work%c, work%scratch, z_zero, work%e)
225 : ! g_surf <- B_{n+1} = B_n A_n^-1 B_n
226 : ! keep B_n, as we need it to compute the matrix product B_n A_n^-1 C_n
227 37416 : CALL parallel_gemm('N', 'N', nrows, nrows, nrows, z_one, work%b, work%scratch, z_zero, g_surf)
228 :
229 : ! scratch <- A_n^-1 * C_n
230 37416 : CALL parallel_gemm('N', 'N', nrows, nrows, nrows, z_one, work%a_inv, work%c, z_zero, work%scratch)
231 : ! D_n = - B_n A_n^-1 C_n
232 37416 : CALL cp_cfm_gemm('N', 'N', nrows, nrows, nrows, z_mone, work%b, work%scratch, z_zero, work%d)
233 : ! we do not need B_n any longer, so the matrix B now holds the B_{n+1} matrix
234 37416 : CALL cp_cfm_to_cfm(g_surf, work%b)
235 : ! C_{n+1} = C_n A_n^-1 C_n
236 37416 : CALL parallel_gemm('N', 'N', nrows, nrows, nrows, z_one, work%c, work%scratch, z_zero, g_surf)
237 37416 : CALL cp_cfm_to_cfm(g_surf, work%c)
238 :
239 : ! A0_{n+1} = A0_n + D_n = A0_n - B_n A_n^-1 C_n
240 37416 : CALL cp_cfm_scale_and_add(z_one, work%a0, z_one, work%d)
241 :
242 : ! A_{n+1} = A0_n + D_n + E_n = A_n - B_n A_n^-1 C_n - C_n A_n^-1 B_n
243 37416 : CALL cp_cfm_scale_and_add(z_one, work%a, z_one, work%d)
244 44026 : CALL cp_cfm_scale_and_add(z_one, work%a, z_one, work%e)
245 : END DO
246 :
247 : ! g_surf = A0_n^-1
248 6610 : CALL cp_cfm_to_cfm(work%a0, g_surf)
249 6610 : CALL cp_cfm_lu_invert(g_surf)
250 :
251 6610 : CALL timestop(handle)
252 6610 : END SUBROUTINE do_sancho
253 :
254 : ! **************************************************************************************************
255 : !> \brief Compute the contact self energy at point 'omega' as
256 : !> self_energy_C = [omega * S_SC0 - KS_SC0] * g_surf_c(omega - v_C) * [omega * S_SC0 - KS_SC0]^T,
257 : !> where C stands for the left (L) / right (R) contact.
258 : !> \param self_energy_c contact self energy (initialised on exit)
259 : !> \param omega energy point where the contact self energy needs to be computed
260 : !> \param g_surf_c contact surface Green's function
261 : !> \param h_sc0 scattering region -- contact off-diagonal block of the Kohn-Sham matrix
262 : !> \param s_sc0 scattering region -- contact off-diagonal block of the overlap matrix
263 : !> \param zwork1 complex work matrix of the same shape as s_sc0
264 : !> \param zwork2 another complex work matrix of the same shape as s_sc0
265 : !> \param transp flag which indicates that transposed matrices (KS_C0S and S_C0S)
266 : !> were actually passed
267 : !> \author Sergey Chulkov
268 : ! **************************************************************************************************
269 11534 : SUBROUTINE negf_contact_self_energy(self_energy_c, omega, g_surf_c, h_sc0, s_sc0, zwork1, zwork2, transp)
270 : TYPE(cp_cfm_type), INTENT(IN) :: self_energy_c
271 : COMPLEX(kind=dp), INTENT(in) :: omega
272 : TYPE(cp_cfm_type), INTENT(IN) :: g_surf_c
273 : TYPE(cp_fm_type), INTENT(IN) :: h_sc0, s_sc0
274 : TYPE(cp_cfm_type), INTENT(IN) :: zwork1, zwork2
275 : LOGICAL, INTENT(in) :: transp
276 :
277 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_contact_self_energy'
278 :
279 : INTEGER :: handle, nao_contact, nao_scattering
280 :
281 11534 : CALL timeset(routineN, handle)
282 :
283 : ! zwork1 = omega * S_SC0 - KS_SC0 if transp = .FALSE., or
284 : ! = omega * S_SC0^T - KS_SC0^T if transp = .TRUE.
285 11534 : CALL cp_fm_to_cfm(msourcer=s_sc0, mtarget=zwork1)
286 11534 : CALL cp_cfm_scale_and_add_fm(omega, zwork1, z_mone, h_sc0)
287 :
288 11534 : IF (transp) THEN
289 1260 : CALL cp_fm_get_info(s_sc0, nrow_global=nao_contact, ncol_global=nao_scattering)
290 :
291 : ! zwork2 = g_surf_c * [omega * S_SC0^T - KS_SC0^T]
292 1260 : CALL parallel_gemm('N', 'N', nao_contact, nao_scattering, nao_contact, z_one, g_surf_c, zwork1, z_zero, zwork2)
293 : ! [omega * S_SC0^T - KS_SC0^T]^T * g_surf_c * [omega * S_SC0^T - KS_SC0^T]
294 1260 : CALL parallel_gemm('T', 'N', nao_scattering, nao_scattering, nao_contact, z_one, zwork1, zwork2, z_zero, self_energy_c)
295 : ELSE
296 10274 : CALL cp_fm_get_info(s_sc0, nrow_global=nao_scattering, ncol_global=nao_contact)
297 :
298 : ! zwork2 = [omega * S_SC0 - KS_SC0] * g_surf_c
299 10274 : CALL parallel_gemm('N', 'N', nao_scattering, nao_contact, nao_contact, z_one, zwork1, g_surf_c, z_zero, zwork2)
300 : ! [omega * S_SC0 - KS_SC0] * g_surf_c * [omega * S_SC0 - KS_SC0]^T
301 10274 : CALL parallel_gemm('N', 'T', nao_scattering, nao_scattering, nao_contact, z_one, zwork2, zwork1, z_zero, self_energy_c)
302 : END IF
303 :
304 11534 : CALL timestop(handle)
305 11534 : END SUBROUTINE negf_contact_self_energy
306 :
307 : ! **************************************************************************************************
308 : !> \brief Compute contact broadening matrix as
309 : !> gamma_C = i (self_energy_c^{ret.} - (self_energy_c^{ret.})^C)
310 : !> \param gamma_c broadening matrix (initialised on exit)
311 : !> \param self_energy_c retarded contact self-energy
312 : !> \author Sergey Chulkov
313 : ! **************************************************************************************************
314 1230 : SUBROUTINE negf_contact_broadening_matrix(gamma_c, self_energy_c)
315 : TYPE(cp_cfm_type), INTENT(IN) :: gamma_c, self_energy_c
316 :
317 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_contact_broadening_matrix'
318 :
319 : INTEGER :: handle
320 :
321 1230 : CALL timeset(routineN, handle)
322 :
323 : ! gamma_contact = i (self_energy_contact^{ret.} - self_energy_contact^{adv.}) .
324 : !
325 : ! With no k-points, the Hamiltonian matrix is real-values, so
326 : ! self_energy_contact^{adv.} = self_energy_contact^{ret.}^C ,
327 : ! The above identity allows us to use a simplified expression for the broadening matrix:
328 : ! gamma_contact = i [self_energy_contact - self_energy_contact^C] .
329 :
330 1230 : CALL cp_cfm_transpose(self_energy_c, 'C', gamma_c)
331 1230 : CALL cp_cfm_scale_and_add(z_mone, gamma_c, z_one, self_energy_c)
332 1230 : CALL cp_cfm_scale(gaussi, gamma_c)
333 :
334 1230 : CALL timestop(handle)
335 1230 : END SUBROUTINE negf_contact_broadening_matrix
336 :
337 : ! **************************************************************************************************
338 : !> \brief Compute the retarded Green's function at point 'omega' as
339 : !> G_S^{ret.} = [ omega * S_S - KS_S - \sum_{contact} self_energy_{contact}^{ret.}]^{-1}.
340 : !> \param g_ret_s complex matrix to store the computed retarded Green's function
341 : !> \param omega energy point where the retarded Green's function needs to be computed
342 : !> \param self_energy_ret_sum sum of contact self-energies
343 : !> \param h_s Kohn-Sham matrix block of the scattering region
344 : !> \param s_s overlap matrix block of the scattering region
345 : !> \param v_hartree_s contribution to the Kohn-Sham matrix from the external potential
346 : !> \author Sergey Chulkov
347 : ! **************************************************************************************************
348 5767 : SUBROUTINE negf_retarded_green_function(g_ret_s, omega, self_energy_ret_sum, h_s, s_s, v_hartree_s)
349 : TYPE(cp_cfm_type), INTENT(IN) :: g_ret_s
350 : COMPLEX(kind=dp), INTENT(in) :: omega
351 : TYPE(cp_cfm_type), INTENT(IN) :: self_energy_ret_sum
352 : TYPE(cp_fm_type), INTENT(IN) :: h_s, s_s
353 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: v_hartree_s
354 :
355 : CHARACTER(len=*), PARAMETER :: routineN = 'negf_retarded_green_function'
356 :
357 : INTEGER :: handle
358 :
359 5767 : CALL timeset(routineN, handle)
360 :
361 : ! g_ret_s = [ omega * S_S - H_S - self_energy_left - self_energy_right]^{-1}
362 : !
363 : ! omega * S_S - H_S - V_Hartree
364 5767 : CALL cp_fm_to_cfm(msourcer=s_s, mtarget=g_ret_s)
365 5767 : CALL cp_cfm_scale_and_add_fm(omega, g_ret_s, z_mone, h_s)
366 5767 : IF (PRESENT(v_hartree_s)) &
367 1597 : CALL cp_cfm_scale_and_add_fm(z_one, g_ret_s, z_one, v_hartree_s)
368 :
369 : ! g_ret_s = [omega * S_S - H_S - \sum_{contact} self_energy_{contact}^{ret.} ]^-1
370 5767 : CALL cp_cfm_scale_and_add(z_one, g_ret_s, z_mone, self_energy_ret_sum)
371 :
372 5767 : CALL cp_cfm_lu_invert(g_ret_s)
373 :
374 5767 : CALL timestop(handle)
375 5767 : END SUBROUTINE negf_retarded_green_function
376 0 : END MODULE negf_green_methods
377 :
|
