Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Static COHSEX band structure in the tensor-based large-cell Gamma GW code:
10 : !> COHSEX DIAGONAL (diagonal energies) or COHSEX QP_UPDATE (canonical energies and MOs)
11 : !> \author Ritaj Tyagi
12 : !> \date Jan 2025
13 : ! **************************************************************************************************
14 : MODULE gw_tensor_large_cell_Gamma_static_cohsex
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE bibliography, ONLY: Hedin1965,&
17 : Hybertsen1986,&
18 : Tyagi2024,&
19 : Tyagi2026,&
20 : cite_reference
21 : USE cell_types, ONLY: cell_type
22 : USE cp_cfm_diag, ONLY: cp_cfm_geeig
23 : USE cp_cfm_types, ONLY: cp_cfm_create,&
24 : cp_cfm_release,&
25 : cp_cfm_to_fm,&
26 : cp_cfm_type
27 : USE cp_control_types, ONLY: dft_control_type
28 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
29 : dbcsr_p_type,&
30 : dbcsr_set
31 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 : dbcsr_deallocate_matrix_set
33 : USE cp_files, ONLY: close_file,&
34 : open_file
35 : USE cp_fm_basic_linalg, ONLY: cp_fm_scale_and_add,&
36 : cp_fm_transpose
37 : USE cp_fm_diag, ONLY: cp_fm_syevx
38 : USE cp_fm_types, ONLY: cp_fm_create,&
39 : cp_fm_get_diag,&
40 : cp_fm_get_info,&
41 : cp_fm_release,&
42 : cp_fm_set_all,&
43 : cp_fm_to_fm,&
44 : cp_fm_type,&
45 : cp_fm_write_formatted
46 : USE cp_log_handling, ONLY: cp_get_default_logger,&
47 : cp_logger_type
48 : USE cp_output_handling, ONLY: cp_p_file,&
49 : cp_print_key_finished_output,&
50 : cp_print_key_should_output,&
51 : cp_print_key_unit_nr
52 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
53 : USE dbt_api, ONLY: dbt_type
54 : USE gw_tensor_large_cell_gamma, ONLY: &
55 : apply_extrapol_factor, compute_MinvVsqrt_Vsqrt, compute_Sigma_c_occ_vir, &
56 : compute_V_k_by_lattice_sum, compute_fm_W_MIC_freq_j_from_chi, create_mat_for_Sigma_c, &
57 : delete_unnecessary_files, destroy_mat_Sigma_c, fm_read, fm_write, get_Sigma_x, &
58 : get_mat_chi_Gamma_tau, print_skipping
59 : USE gw_utils, ONLY: de_init_bs_env
60 : USE gw_utils_fm, ONLY: cfm_contract_ABA,&
61 : fm_contract_ABA
62 : USE input_constants, ONLY: cohsex_qp_update
63 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
64 : section_vals_type,&
65 : section_vals_val_get
66 : USE kinds, ONLY: default_path_length,&
67 : default_string_length,&
68 : dp
69 : USE machine, ONLY: m_walltime
70 : USE parallel_gemm_api, ONLY: parallel_gemm
71 : USE particle_list_types, ONLY: particle_list_type
72 : USE particle_types, ONLY: particle_type
73 : USE physcon, ONLY: evolt
74 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
75 : USE post_scf_bandstructure_utils, ONLY: cfm_ikp_from_fm_Gamma,&
76 : get_all_VBM_CBM_bandgaps
77 : USE pw_env_types, ONLY: pw_env_get,&
78 : pw_env_type
79 : USE pw_pool_types, ONLY: pw_pool_type
80 : USE pw_types, ONLY: pw_c1d_gs_type,&
81 : pw_r3d_rs_type
82 : USE qs_collocate_density, ONLY: calculate_wavefunction
83 : USE qs_environment_types, ONLY: get_qs_env,&
84 : qs_environment_type
85 : USE qs_kind_types, ONLY: qs_kind_type
86 : USE qs_mo_io, ONLY: write_mo_set_to_output_unit
87 : USE qs_mo_types, ONLY: allocate_mo_set,&
88 : deallocate_mo_set,&
89 : init_mo_set,&
90 : mo_set_type
91 : USE qs_subsys_types, ONLY: qs_subsys_get,&
92 : qs_subsys_type
93 : #include "./base/base_uses.f90"
94 :
95 : IMPLICIT NONE
96 :
97 : PRIVATE
98 :
99 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_tensor_large_cell_Gamma_static_cohsex'
100 :
101 : PUBLIC :: gw_calc_tensor_large_cell_Gamma_static_cohsex
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief Perform static COHSEX band structure calculation
107 : !> \param qs_env ...
108 : !> \param bs_env ...
109 : !> \par History
110 : !> * 12.2024 created [Ritaj Tyagi]
111 : ! **************************************************************************************************
112 8 : SUBROUTINE gw_calc_tensor_large_cell_Gamma_static_cohsex(qs_env, bs_env)
113 : TYPE(qs_environment_type), POINTER :: qs_env
114 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
115 :
116 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gw_calc_tensor_large_cell_Gamma_static_cohsex'
117 :
118 : INTEGER :: handle
119 8 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_c_Gamma, fm_Sigma_x_Gamma, &
120 8 : fm_W_MIC_freq_zero
121 :
122 8 : CALL timeset(routineN, handle)
123 :
124 : ! static COHSEX; with QP_UPDATE the quasiparticle update; scaled self-energy if requested
125 8 : CALL cite_reference(Hedin1965)
126 8 : CALL cite_reference(Hybertsen1986)
127 8 : IF (bs_env%cohsex_flavor == cohsex_qp_update) CALL cite_reference(Tyagi2024)
128 8 : IF (bs_env%cohsex_scale_occ /= 1.0_dp .OR. bs_env%cohsex_scale_virt /= 1.0_dp) THEN
129 8 : CALL cite_reference(Tyagi2026)
130 : END IF
131 :
132 8 : IF (bs_env%Sigma_c_freq_zero_exists) THEN
133 :
134 : ! restart: Σ^c_λσ(iω=0,k=0) from file; χ and W are not needed
135 0 : CALL read_Sigma_c_freq_zero(bs_env, fm_Sigma_c_Gamma)
136 :
137 : ELSE
138 :
139 : ! G^occ_µλ(i|τ|,k=0) = sum_n^occ C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
140 : ! G^vir_µλ(i|τ|,k=0) = sum_n^vir C_µn(k=0) e^(-|(ϵ_nk=0-ϵ_F)τ|) C_λn(k=0)
141 : ! χ_PQ(iτ,k=0) = sum_λν [sum_µ (µν|P) G^occ_µλ(i|τ|)] [sum_σ (σλ|Q) G^vir_σν(i|τ|)]
142 8 : CALL get_mat_chi_Gamma_tau(bs_env, qs_env, bs_env%mat_chi_Gamma_tau)
143 :
144 : ! χ_PQ(iτ,k=0) -> χ_PQ(iω=0,k=0) -> χ_PQ(iω=0,k) -> ε_PQ(iω=0,k) -> W_PQ(iω=0,k)
145 : ! -> W^MIC_PQ(iω=0) -> M^-1*W^MIC(iω=0)*M^-1
146 8 : CALL get_W_MIC_freq_zero(bs_env, qs_env, bs_env%mat_chi_Gamma_tau, fm_W_MIC_freq_zero)
147 :
148 : ! Σ^c_λσ(iω=0,k=0) = 1/2 sum_n sgn(ϵ_n - ϵ_F) sum_νQ [sum_P (νσ|P) W^MIC_PQ(iω=0)]
149 : ! [sum_µ (λµ|Q) C_µn C_νn]; written to the Σ^c(iω=0) restart file
150 8 : CALL get_Sigma_c(bs_env, qs_env, fm_W_MIC_freq_zero, fm_Sigma_c_Gamma)
151 :
152 : END IF
153 :
154 : ! D_µν = sum_n^occ C_µn(k=0) C_νn(k=0), V^trunc_PQ = sum_cell_R <phi_P,0|V^trunc|phi_Q,R>
155 : ! Σ^x_λσ(k=0) = sum_νQ [sum_P (νσ|P) V^trunc_PQ] [sum_µ (λµ|Q) D_µν)]
156 8 : CALL get_Sigma_x(bs_env, qs_env, fm_Sigma_x_Gamma)
157 :
158 : ! Σ^c_λσ(k=0) -> Σ^c_nn(ϵ,k); ϵ_nk^GW = ϵ_nk^DFT + Σ^c_nn(ϵ,k) + Σ^x_nn(k) - v^xc_nn(k)
159 8 : CALL compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma)
160 :
161 8 : CALL de_init_bs_env(qs_env, bs_env)
162 :
163 8 : CALL timestop(handle)
164 :
165 8 : END SUBROUTINE gw_calc_tensor_large_cell_Gamma_static_cohsex
166 :
167 : ! **************************************************************************************************
168 : !> \brief Static screened Coulomb interaction at iω = 0:
169 : !> χ_PQ(iω=0,k=0) = sum_i w^0_i χ_PQ(iτ_i,k=0) (w^0: time weights at zero frequency)
170 : !> χ_PQ(iω=0,k=0) -> χ_PQ(iω=0,k) -> ε_PQ(iω=0,k) -> W_PQ(iω=0,k) -> W^MIC_PQ(iω=0)
171 : !> as for one frequency point of G0W0 (compute_fm_W_MIC_freq_j_from_chi); then the
172 : !> k-point extrapolation factor and M^-1(k=0) W^MIC(iω=0) M^-1(k=0)
173 : !> \param bs_env ...
174 : !> \param qs_env ...
175 : !> \param mat_chi_Gamma_tau ...
176 : !> \param fm_W_MIC_freq_zero ...
177 : ! **************************************************************************************************
178 8 : SUBROUTINE get_W_MIC_freq_zero(bs_env, qs_env, mat_chi_Gamma_tau, fm_W_MIC_freq_zero)
179 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
180 : TYPE(qs_environment_type), POINTER :: qs_env
181 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mat_chi_Gamma_tau
182 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_freq_zero
183 :
184 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_W_MIC_freq_zero'
185 :
186 : INTEGER :: handle, i_t, ikp, ikp_batch, ikp_in_batch
187 : REAL(KIND=dp) :: t1
188 : TYPE(cp_cfm_type) :: cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp
189 8 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:) :: cfm_V_kp
190 :
191 8 : CALL timeset(routineN, handle)
192 :
193 : ! χ_PQ(iω=0,k=0) = sum_i w^0_i χ_PQ(iτ_i,k=0) -> bs_env%fm_chi_Gamma_freq
194 8 : CALL dbcsr_set(bs_env%mat_RI_RI%matrix, 0.0_dp)
195 60 : DO i_t = 1, bs_env%num_time_freq_points
196 : CALL dbcsr_add(bs_env%mat_RI_RI%matrix, mat_chi_Gamma_tau(i_t)%matrix, 1.0_dp, &
197 60 : bs_env%time_frequency_grid%time_weights_at_zero_frequency(i_t))
198 : END DO
199 8 : CALL copy_dbcsr_to_fm(bs_env%mat_RI_RI%matrix, bs_env%fm_chi_Gamma_freq)
200 8 : CALL dbcsr_deallocate_matrix_set(mat_chi_Gamma_tau)
201 :
202 16 : ALLOCATE (fm_W_MIC_freq_zero(1))
203 8 : CALL cp_fm_create(fm_W_MIC_freq_zero(1), bs_env%fm_RI_RI%matrix_struct)
204 8 : CALL cp_fm_set_all(fm_W_MIC_freq_zero(1), 0.0_dp)
205 :
206 62 : DO ikp_batch = 1, bs_env%num_chi_eps_W_batches
207 :
208 54 : t1 = m_walltime()
209 :
210 : ! V_PQ(k) = sum_R e^(ikR) <phi_P, cell 0 | 1/r | phi_Q, cell R>
211 54 : CALL compute_V_k_by_lattice_sum(bs_env, qs_env, cfm_V_kp, ikp_batch)
212 :
213 270 : DO ikp_in_batch = 1, bs_env%nkp_chi_eps_W_batch
214 :
215 216 : ikp = (ikp_batch - 1)*bs_env%nkp_chi_eps_W_batch + ikp_in_batch
216 :
217 216 : IF (ikp > bs_env%nkp_chi_eps_W_orig_plus_extra) CYCLE
218 :
219 : CALL compute_MinvVsqrt_Vsqrt(bs_env, qs_env, cfm_V_kp, &
220 206 : cfm_V_sqrt_ikp, cfm_M_inv_V_sqrt_ikp, ikp)
221 :
222 206 : CALL bs_env%para_env%sync()
223 206 : CALL cp_cfm_release(cfm_V_kp(ikp))
224 :
225 : ! W^MIC contribution of k-point ikp at iω = 0 (j_w = 1: first frequency point for the
226 : ! approximate k-point extrapolation)
227 : CALL compute_fm_W_MIC_freq_j_from_chi(bs_env, qs_env, bs_env%fm_W_MIC_freq, 1, ikp, &
228 206 : cfm_M_inv_V_sqrt_ikp, cfm_V_sqrt_ikp)
229 :
230 270 : CALL cp_fm_scale_and_add(1.0_dp, fm_W_MIC_freq_zero(1), 1.0_dp, bs_env%fm_W_MIC_freq)
231 :
232 : END DO ! ikp_in_batch
233 :
234 54 : DEALLOCATE (cfm_V_kp)
235 :
236 62 : IF (bs_env%unit_nr > 0) THEN
237 : WRITE (bs_env%unit_nr, '(T2,A,I10,A,I3,A,F10.1,A)') &
238 27 : 'Computed W(iω=0,k) for k-point batch', &
239 27 : ikp_batch, ' /', bs_env%num_chi_eps_W_batches, &
240 54 : ', Execution time', m_walltime() - t1, ' s'
241 : END IF
242 :
243 : END DO ! ikp_batch
244 :
245 8 : IF (bs_env%approx_kp_extrapol) THEN
246 0 : CALL apply_extrapol_factor(bs_env, fm_W_MIC_freq_zero)
247 : END IF
248 :
249 : ! M^-1(k=0) W^MIC(iω=0) M^-1(k=0)
250 8 : CALL fm_contract_ABA(bs_env%fm_Minv_Gamma, fm_W_MIC_freq_zero(1))
251 :
252 8 : CALL cp_cfm_release(cfm_M_inv_V_sqrt_ikp)
253 8 : CALL cp_cfm_release(cfm_V_sqrt_ikp)
254 :
255 8 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
256 :
257 8 : CALL timestop(handle)
258 :
259 16 : END SUBROUTINE get_W_MIC_freq_zero
260 :
261 : ! **************************************************************************************************
262 : !> \brief Static COHSEX correlation self-energy at the Γ point in the AO basis,
263 : !> Σ^c_λσ = 1/2 sum_n sgn(ϵ_n - ϵ_F) sum_νQ [sum_P (νσ|P) W^MIC_PQ(iω=0)] [sum_µ (λµ|Q) C_µn C_νn]
264 : !> = 1/2 (Σ^neg_λσ(iτ=0) + Σ^pos_λσ(iτ=0)) with W = W^MIC(iω=0),
265 : !> i.e. the occupied and virtual self-energy parts of G0W0 at τ = 0 (compute_Sigma_c_occ_vir),
266 : !> symmetrized
267 : !> \param bs_env ...
268 : !> \param qs_env ...
269 : !> \param fm_W_MIC_freq_zero ...
270 : !> \param fm_Sigma_c_Gamma ...
271 : ! **************************************************************************************************
272 8 : SUBROUTINE get_Sigma_c(bs_env, qs_env, fm_W_MIC_freq_zero, fm_Sigma_c_Gamma)
273 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
274 : TYPE(qs_environment_type), POINTER :: qs_env
275 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_W_MIC_freq_zero, fm_Sigma_c_Gamma
276 :
277 : CHARACTER(LEN=*), PARAMETER :: routineN = 'get_Sigma_c'
278 :
279 : INTEGER :: handle, ispin
280 : REAL(KIND=dp) :: t1
281 8 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: mat_Sigma_neg, mat_Sigma_pos
282 200 : TYPE(dbt_type) :: t_2c_Gocc, t_2c_Gvir, t_2c_Sigma_neg, &
283 200 : t_2c_Sigma_pos, t_2c_W, t_3c_x_W
284 :
285 8 : CALL timeset(routineN, handle)
286 :
287 : CALL create_mat_for_Sigma_c(bs_env, t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg, &
288 8 : t_2c_Sigma_pos, t_3c_x_W, mat_Sigma_neg, mat_Sigma_pos)
289 :
290 32 : ALLOCATE (fm_Sigma_c_Gamma(bs_env%n_spin))
291 16 : DO ispin = 1, bs_env%n_spin
292 :
293 8 : t1 = m_walltime()
294 :
295 : ! Σ^neg(iτ=0) = -G^occ(0) W^MIC(iω=0), Σ^pos(iτ=0) = G^vir(0) W^MIC(iω=0)
296 : CALL compute_Sigma_c_occ_vir(bs_env, qs_env, 0.0_dp, ispin, fm_W_MIC_freq_zero(1), &
297 : t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg, &
298 : t_2c_Sigma_pos, t_3c_x_W, mat_Sigma_neg(1, ispin)%matrix, &
299 8 : mat_Sigma_pos(1, ispin)%matrix, count_skipped=.FALSE.)
300 :
301 : ! Σ^c = 1/2 (Σ^neg + Σ^pos), symmetrized
302 8 : CALL cp_fm_create(fm_Sigma_c_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
303 8 : CALL copy_dbcsr_to_fm(mat_Sigma_neg(1, ispin)%matrix, fm_Sigma_c_Gamma(ispin))
304 8 : CALL copy_dbcsr_to_fm(mat_Sigma_pos(1, ispin)%matrix, bs_env%fm_work_mo(1))
305 8 : CALL cp_fm_scale_and_add(0.5_dp, fm_Sigma_c_Gamma(ispin), 0.5_dp, bs_env%fm_work_mo(1))
306 8 : CALL cp_fm_transpose(fm_Sigma_c_Gamma(ispin), bs_env%fm_work_mo(1))
307 8 : CALL cp_fm_scale_and_add(0.5_dp, fm_Sigma_c_Gamma(ispin), 0.5_dp, bs_env%fm_work_mo(1))
308 :
309 : ! restart file of Σ^c(iω=0,k=0), one per spin
310 8 : CALL fm_write(fm_Sigma_c_Gamma(ispin), ispin, bs_env%Sigma_c_freq_zero_name, qs_env)
311 :
312 16 : IF (bs_env%unit_nr > 0) THEN
313 : WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
314 4 : 'Computed Σ^c(k=0),', ' Execution time', m_walltime() - t1, ' s'
315 : END IF
316 :
317 : END DO
318 :
319 8 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
320 :
321 8 : CALL print_skipping(bs_env)
322 :
323 : CALL destroy_mat_Sigma_c(t_2c_Gocc, t_2c_Gvir, t_2c_W, t_2c_Sigma_neg, t_2c_Sigma_pos, &
324 8 : t_3c_x_W, fm_W_MIC_freq_zero, mat_Sigma_neg, mat_Sigma_pos)
325 :
326 8 : CALL delete_unnecessary_files(bs_env)
327 :
328 8 : CALL timestop(handle)
329 :
330 24 : END SUBROUTINE get_Sigma_c
331 :
332 : ! **************************************************************************************************
333 : !> \brief Read the static COHSEX restart file Σ^c_λσ(iω=0,k=0), one per spin
334 : !> \param bs_env ...
335 : !> \param fm_Sigma_c_Gamma Σ^c_λσ(iω=0,k=0) per spin, created here
336 : ! **************************************************************************************************
337 0 : SUBROUTINE read_Sigma_c_freq_zero(bs_env, fm_Sigma_c_Gamma)
338 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
339 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_c_Gamma
340 :
341 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_Sigma_c_freq_zero'
342 :
343 : INTEGER :: handle, ispin
344 : REAL(KIND=dp) :: t1
345 :
346 0 : CALL timeset(routineN, handle)
347 :
348 0 : ALLOCATE (fm_Sigma_c_Gamma(bs_env%n_spin))
349 0 : DO ispin = 1, bs_env%n_spin
350 0 : t1 = m_walltime()
351 0 : CALL cp_fm_create(fm_Sigma_c_Gamma(ispin), bs_env%fm_s_Gamma%matrix_struct)
352 0 : CALL fm_read(fm_Sigma_c_Gamma(ispin), bs_env, bs_env%Sigma_c_freq_zero_name, ispin)
353 0 : IF (bs_env%unit_nr > 0) THEN
354 : WRITE (bs_env%unit_nr, '(T2,A,T55,A,F10.1,A)') &
355 0 : 'Read Σ^c(k=0) from file,', ' Execution time', m_walltime() - t1, ' s'
356 : END IF
357 : END DO
358 :
359 0 : IF (bs_env%unit_nr > 0) WRITE (bs_env%unit_nr, '(A)') ' '
360 :
361 0 : CALL timestop(handle)
362 :
363 0 : END SUBROUTINE read_Sigma_c_freq_zero
364 :
365 : ! **************************************************************************************************
366 : !> \brief ...
367 : !> \param bs_env ...
368 : !> \param qs_env ...
369 : !> \param fm_Sigma_x_Gamma ...
370 : !> \param fm_Sigma_c_Gamma ...
371 : ! **************************************************************************************************
372 48 : SUBROUTINE compute_QP_energies(bs_env, qs_env, fm_Sigma_x_Gamma, fm_Sigma_c_Gamma)
373 :
374 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
375 : TYPE(qs_environment_type), POINTER :: qs_env
376 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: fm_Sigma_x_Gamma, fm_Sigma_c_Gamma
377 :
378 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_QP_energies'
379 :
380 : INTEGER :: handle, ikp, ispin
381 : TYPE(cp_cfm_type) :: cfm_ks_ikp, cfm_mos_ikp, cfm_s_ikp, &
382 : cfm_work_ikp
383 : TYPE(cp_fm_type) :: fm_Sigma_c_ikp_mo, fm_Sigma_x_ikp_mo, &
384 : fm_V_xc_ikp_mo
385 :
386 8 : CALL timeset(routineN, handle)
387 :
388 8 : CALL cp_cfm_create(cfm_mos_ikp, bs_env%fm_s_Gamma%matrix_struct)
389 8 : CALL cp_cfm_create(cfm_work_ikp, bs_env%fm_s_Gamma%matrix_struct)
390 :
391 : ! MO-basis operator matrices at k_i (n_ao x n_ao), reused for every k-point
392 8 : CALL cp_fm_create(fm_V_xc_ikp_mo, bs_env%fm_s_Gamma%matrix_struct)
393 8 : CALL cp_fm_create(fm_Sigma_x_ikp_mo, bs_env%fm_s_Gamma%matrix_struct)
394 8 : CALL cp_fm_create(fm_Sigma_c_ikp_mo, bs_env%fm_s_Gamma%matrix_struct)
395 :
396 16 : DO ispin = 1, bs_env%n_spin
397 :
398 34 : DO ikp = 1, bs_env%nkp_bs_and_DOS
399 :
400 : ! 1. get H^KS_µν(k_i) from H^KS_µν(k=0)
401 : CALL cfm_ikp_from_fm_Gamma(cfm_ks_ikp, bs_env%fm_ks_Gamma(ispin), &
402 18 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
403 :
404 : ! 2. get S_µν(k_i) from S_µν(k=0)
405 : CALL cfm_ikp_from_fm_Gamma(cfm_s_ikp, bs_env%fm_s_Gamma, &
406 18 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
407 :
408 : ! 3. Diagonalize (Roothaan-Hall): H_KS(k_i)*C(k_i) = S(k_i)*C(k_i)*ϵ(k_i)
409 : CALL cp_cfm_geeig(cfm_ks_ikp, cfm_s_ikp, cfm_mos_ikp, &
410 18 : bs_env%eigenval_scf(:, ikp, ispin), cfm_work_ikp)
411 :
412 : ! Reuse the KS buffer for the full MO operators. Retain their real parts
413 : ! at the existing real COHSEX Hamiltonian boundary.
414 : ! 4. V^xc_µν(k=0) -> V^xc_µν(k_i) -> V^xc_mn(k_i)
415 : CALL cfm_ikp_from_fm_Gamma(cfm_ks_ikp, bs_env%fm_V_xc_Gamma(ispin), &
416 18 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
417 18 : CALL cfm_contract_ABA(cfm_mos_ikp, cfm_ks_ikp, workspace=cfm_work_ikp)
418 18 : CALL cp_cfm_to_fm(cfm_ks_ikp, fm_V_xc_ikp_mo)
419 :
420 : ! 5. Σ^x_µν(k=0) -> Σ^x_µν(k_i) -> Σ^x_mn(k_i)
421 : CALL cfm_ikp_from_fm_Gamma(cfm_ks_ikp, fm_Sigma_x_Gamma(ispin), &
422 18 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
423 18 : CALL cfm_contract_ABA(cfm_mos_ikp, cfm_ks_ikp, workspace=cfm_work_ikp)
424 18 : CALL cp_cfm_to_fm(cfm_ks_ikp, fm_Sigma_x_ikp_mo)
425 :
426 : ! 6. Σ^c_µν(k=0) -> Σ^c_µν(k_i) -> Σ^c_mn(k_i)
427 : CALL cfm_ikp_from_fm_Gamma(cfm_ks_ikp, fm_Sigma_c_Gamma(ispin), &
428 18 : ikp, qs_env, bs_env%kpoints_DOS, "ORB")
429 18 : CALL cfm_contract_ABA(cfm_mos_ikp, cfm_ks_ikp, workspace=cfm_work_ikp)
430 18 : CALL cp_cfm_to_fm(cfm_ks_ikp, fm_Sigma_c_ikp_mo)
431 :
432 : ! 7. build static COHSEX QP Hamiltonian, diagonalize, and print
433 : CALL diagonalize_and_print(bs_env, qs_env, fm_Sigma_x_ikp_mo, fm_V_xc_ikp_mo, &
434 26 : fm_Sigma_c_ikp_mo, ikp, ispin)
435 :
436 : END DO ! ikp_DOS
437 :
438 : END DO ! ispin
439 :
440 8 : CALL get_all_VBM_CBM_bandgaps(bs_env)
441 :
442 8 : CALL cp_fm_release(fm_Sigma_x_Gamma)
443 8 : CALL cp_fm_release(fm_Sigma_c_Gamma)
444 8 : CALL cp_fm_release(fm_V_xc_ikp_mo)
445 8 : CALL cp_fm_release(fm_Sigma_x_ikp_mo)
446 8 : CALL cp_fm_release(fm_Sigma_c_ikp_mo)
447 8 : CALL cp_cfm_release(cfm_ks_ikp)
448 8 : CALL cp_cfm_release(cfm_s_ikp)
449 8 : CALL cp_cfm_release(cfm_mos_ikp)
450 8 : CALL cp_cfm_release(cfm_work_ikp)
451 :
452 8 : CALL timestop(handle)
453 :
454 8 : END SUBROUTINE compute_QP_energies
455 :
456 : ! **************************************************************************************************
457 : !> \brief Build the static COHSEX quasiparticle Hamiltonian H^QP at k-point ikp and print the
458 : !> diagonal energies ϵ^DG = H^QP_nn (COHSEX DIAGONAL); with COHSEX QP_UPDATE also
459 : !> diagonalize H^QP and print the canonical energies ϵ^CN and the canonical MOs
460 : !> \param bs_env ...
461 : !> \param qs_env ...
462 : !> \param fm_Sigma_x_ikp_mo ...
463 : !> \param fm_V_xc_ikp_mo ...
464 : !> \param fm_Sigma_c_ikp_mo ...
465 : !> \param ikp ...
466 : !> \param ispin ...
467 : ! **************************************************************************************************
468 18 : SUBROUTINE diagonalize_and_print(bs_env, qs_env, fm_Sigma_x_ikp_mo, fm_V_xc_ikp_mo, &
469 : fm_Sigma_c_ikp_mo, ikp, ispin)
470 :
471 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
472 : TYPE(qs_environment_type), POINTER :: qs_env
473 : TYPE(cp_fm_type) :: fm_Sigma_x_ikp_mo, fm_V_xc_ikp_mo, &
474 : fm_Sigma_c_ikp_mo
475 : INTEGER :: ikp, ispin
476 :
477 : CHARACTER(LEN=*), PARAMETER :: routineN = 'diagonalize_and_print'
478 :
479 : CHARACTER(len=3) :: occ_vir
480 : CHARACTER(len=default_string_length) :: fname
481 : INTEGER :: handle, i_col, i_row, ikp_for_print, &
482 : iunit, j_col, j_row, m, n_mo, &
483 : ncol_local, nrow_local
484 18 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
485 : LOGICAL :: is_bandstruc_kpoint, print_DOS_kpoints, &
486 : print_ikp, qp_update
487 : REAL(KIND=dp) :: val
488 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: diag_Sigma_c, diag_Sigma_x, diag_SP, &
489 : diag_V_xc, SP_eigval
490 : TYPE(cp_fm_type) :: fm_SP_mat, SP_eigvec, SP_eigvec_ao
491 :
492 18 : CALL timeset(routineN, handle)
493 :
494 18 : n_mo = bs_env%n_ao
495 : ! COHSEX QP_UPDATE: diagonalize H^QP (canonical energies and MOs); DIAGONAL: ϵ^DG only
496 18 : qp_update = (bs_env%cohsex_flavor == cohsex_qp_update)
497 :
498 18 : CALL cp_fm_create(fm_SP_mat, bs_env%fm_s_Gamma%matrix_struct)
499 :
500 54 : ALLOCATE (SP_eigval(n_mo))
501 18 : SP_eigval(:) = 0.0_dp
502 :
503 : ! Build the static COHSEX QP Hamiltonian directly in the distributed fm by looping
504 : ! over local elements only (global indices from row_indices/col_indices):
505 : ! occ-occ block : Σ^x - v^xc + α Σ^c (α = COHSEX_SCALE_OCC)
506 : ! vir-vir block : Σ^x - v^xc + α Σ^c (α = COHSEX_SCALE_VIRT)
507 : ! occ-vir / vir-occ : 0, i.e. occupied and virtual states are not mixed
508 : ! diagonal additionally gets + ϵ^DFT_nk
509 : CALL cp_fm_get_info(fm_SP_mat, nrow_local=nrow_local, ncol_local=ncol_local, &
510 18 : row_indices=row_indices, col_indices=col_indices)
511 :
512 : !$OMP PARALLEL DO COLLAPSE(2) DEFAULT(NONE) &
513 : !$OMP SHARED(ncol_local, nrow_local, col_indices, row_indices, bs_env, ispin, ikp, &
514 : !$OMP fm_Sigma_x_ikp_mo, fm_V_xc_ikp_mo, fm_Sigma_c_ikp_mo, fm_SP_mat) &
515 18 : !$OMP PRIVATE(i_col, i_row, j_col, j_row, val)
516 : DO i_col = 1, ncol_local
517 : DO i_row = 1, nrow_local
518 : j_col = col_indices(i_col)
519 : j_row = row_indices(i_row)
520 : IF (j_row <= bs_env%n_occ(ispin) .AND. j_col <= bs_env%n_occ(ispin)) THEN
521 : val = fm_Sigma_x_ikp_mo%local_data(i_row, i_col) &
522 : - fm_V_xc_ikp_mo%local_data(i_row, i_col) &
523 : + bs_env%cohsex_scale_occ*fm_Sigma_c_ikp_mo%local_data(i_row, i_col)
524 : ELSE IF (j_row > bs_env%n_occ(ispin) .AND. j_col > bs_env%n_occ(ispin)) THEN
525 : val = fm_Sigma_x_ikp_mo%local_data(i_row, i_col) &
526 : - fm_V_xc_ikp_mo%local_data(i_row, i_col) &
527 : + bs_env%cohsex_scale_virt*fm_Sigma_c_ikp_mo%local_data(i_row, i_col)
528 : ELSE
529 : val = 0.0_dp
530 : END IF
531 : IF (j_row == j_col) val = val + bs_env%eigenval_scf(j_row, ikp, ispin)
532 : fm_SP_mat%local_data(i_row, i_col) = val
533 : END DO
534 : END DO
535 : !$OMP END PARALLEL DO
536 :
537 : ! diagonal elements (for the per-MO QP printout) before diagonalization
538 90 : ALLOCATE (diag_SP(n_mo), diag_Sigma_c(n_mo), diag_Sigma_x(n_mo), diag_V_xc(n_mo))
539 18 : CALL cp_fm_get_diag(fm_SP_mat, diag_SP)
540 18 : CALL cp_fm_get_diag(fm_Sigma_c_ikp_mo, diag_Sigma_c)
541 18 : CALL cp_fm_get_diag(fm_Sigma_x_ikp_mo, diag_Sigma_x)
542 18 : CALL cp_fm_get_diag(fm_V_xc_ikp_mo, diag_V_xc)
543 :
544 18 : IF (qp_update) THEN
545 : ! canonical static COHSEX: H^QP U = U ϵ^CN (the occ-virt blocks are zero)
546 14 : CALL cp_fm_create(SP_eigvec, bs_env%fm_s_Gamma%matrix_struct)
547 14 : CALL cp_fm_create(SP_eigvec_ao, bs_env%fm_s_Gamma%matrix_struct)
548 14 : CALL cp_fm_syevx(matrix=fm_SP_mat, eigenvectors=SP_eigvec, eigenvalues=SP_eigval)
549 14 : IF (ikp == bs_env%nkp_bs_and_DOS) THEN
550 6 : CALL fm_write_formatted(SP_eigvec, 0, "Singlepole_Eigenvector_Gamma", qs_env)
551 : END IF
552 :
553 : ! transform QP eigenvectors from the KS-MO basis to the AO basis
554 : CALL parallel_gemm(transa="N", transb="N", m=n_mo, n=n_mo, k=n_mo, alpha=1.0_dp, &
555 : matrix_a=bs_env%fm_mo_coeff_Gamma(ispin), matrix_b=SP_eigvec, &
556 14 : beta=0.0_dp, matrix_c=SP_eigvec_ao)
557 : ELSE
558 : ! diagonal static COHSEX: ϵ^DG = H^QP_nn, no diagonalization
559 96 : SP_eigval(:) = diag_SP(:)
560 : END IF
561 :
562 412 : bs_env%eigenval_GW(:, ikp, ispin) = SP_eigval(:)
563 18 : IF (ALLOCATED(bs_env%eigenval_G0W0)) THEN
564 412 : bs_env%eigenval_G0W0(:, ikp, ispin) = SP_eigval(:)
565 : END IF
566 412 : DO m = 1, n_mo
567 : bs_env%eigenval_HF(m, ikp, ispin) = bs_env%eigenval_scf(m, ikp, ispin) + &
568 412 : diag_Sigma_x(m) - diag_V_xc(m)
569 : END DO
570 :
571 : ! only print eigenvalues of DOS k-points in case no bandstructure path has been given;
572 : ! in kpoints_DOS the last nkp_only_bs entries are the bandstructure-path k-points
573 18 : print_DOS_kpoints = (bs_env%nkp_only_bs <= 0)
574 18 : is_bandstruc_kpoint = (ikp > bs_env%nkp_only_DOS)
575 18 : print_ikp = print_DOS_kpoints .OR. is_bandstruc_kpoint
576 :
577 : !! Printing Quasiparticle energies to a file
578 18 : IF (bs_env%para_env%is_source() .AND. print_ikp) THEN
579 :
580 9 : IF (print_DOS_kpoints) THEN
581 9 : ikp_for_print = ikp
582 : ELSE
583 0 : ikp_for_print = ikp - bs_env%nkp_only_DOS
584 : END IF
585 :
586 9 : fname = "Quasiparticle_energies_singlepole"
587 9 : IF (ikp_for_print == 1 .AND. ispin == 1) THEN
588 4 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="REPLACE", file_action="WRITE")
589 : ELSE
590 : CALL open_file(TRIM(fname), unit_number=iunit, file_status="OLD", file_action="WRITE", &
591 5 : file_position="APPEND")
592 : END IF
593 :
594 9 : WRITE (iunit, "(A)") " "
595 9 : WRITE (iunit, "(A10,I7,A25,3F10.4)") "kpoint: ", ikp_for_print, "coordinate: ", &
596 45 : bs_env%kpoints_DOS%xkp(:, ikp)
597 9 : WRITE (iunit, "(A)") " "
598 9 : IF (qp_update) THEN
599 7 : WRITE (iunit, "(A5,A24,2A17,A16,A18,A18)") "n", "ϵ_nk^DFT (eV)", "Σ^c_nk (eV)", &
600 14 : "Σ^x_nk (eV)", "v_nk^xc (eV)", "ϵ_nk^DG (eV)", "ϵ_nk^CN (eV)"
601 : ELSE
602 2 : WRITE (iunit, "(A5,A24,2A17,A16,A18)") "n", "ϵ_nk^DFT (eV)", "Σ^c_nk (eV)", &
603 4 : "Σ^x_nk (eV)", "v_nk^xc (eV)", "ϵ_nk^DG (eV)"
604 : END IF
605 9 : WRITE (iunit, "(A)") " "
606 :
607 206 : DO m = 1, n_mo
608 197 : IF (m <= bs_env%n_occ(ispin)) occ_vir = 'occ'
609 197 : IF (m > bs_env%n_occ(ispin)) occ_vir = 'vir'
610 206 : IF (qp_update) THEN
611 151 : WRITE (iunit, "(I5,3A,4F16.3,F17.3,F17.3)") m, ' (', occ_vir, ') ', &
612 151 : bs_env%eigenval_scf(m, ikp, ispin)*evolt, &
613 151 : diag_Sigma_c(m)*evolt, &
614 151 : diag_Sigma_x(m)*evolt, &
615 151 : diag_V_xc(m)*evolt, &
616 151 : diag_SP(m)*evolt, &
617 302 : SP_eigval(m)*evolt
618 : ELSE
619 46 : WRITE (iunit, "(I5,3A,4F16.3,F17.3)") m, ' (', occ_vir, ') ', &
620 46 : bs_env%eigenval_scf(m, ikp, ispin)*evolt, &
621 46 : diag_Sigma_c(m)*evolt, &
622 46 : diag_Sigma_x(m)*evolt, &
623 46 : diag_V_xc(m)*evolt, &
624 92 : diag_SP(m)*evolt
625 : END IF
626 : END DO
627 :
628 9 : WRITE (iunit, "(A)") " "
629 :
630 9 : CALL close_file(iunit)
631 : END IF
632 :
633 18 : IF (qp_update) THEN
634 : ! canonical MO energies/coefficients and MO cubes (DFT%PRINT%MO, DFT%PRINT%MO_CUBES),
635 : ! Gamma point only
636 14 : IF (ikp == 1) THEN
637 6 : CALL print_static_cohsex_mos(bs_env, qs_env, SP_eigvec_ao, SP_eigval, ispin)
638 : END IF
639 14 : CALL cp_fm_release(SP_eigvec)
640 14 : CALL cp_fm_release(SP_eigvec_ao)
641 : END IF
642 :
643 18 : CALL cp_fm_release(fm_SP_mat)
644 18 : DEALLOCATE (SP_eigval, diag_SP, diag_Sigma_c, diag_Sigma_x, diag_V_xc)
645 :
646 18 : CALL timestop(handle)
647 :
648 36 : END SUBROUTINE diagonalize_and_print
649 :
650 : ! **************************************************************************************************
651 : !> \brief Print static COHSEX QP MO energies/occupations/coefficients and MO cubes at the Gamma point,
652 : !> reusing the standard DFT%PRINT%MO / DFT%PRINT%MO_CUBES print keys.
653 : !> \param bs_env ...
654 : !> \param qs_env ...
655 : !> \param fm_qp_mo_ao ...
656 : !> \param qp_eigval ...
657 : !> \param ispin ...
658 : ! **************************************************************************************************
659 6 : SUBROUTINE print_static_cohsex_mos(bs_env, qs_env, fm_qp_mo_ao, qp_eigval, ispin)
660 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
661 : TYPE(qs_environment_type), POINTER :: qs_env
662 : TYPE(cp_fm_type) :: fm_qp_mo_ao
663 : REAL(KIND=dp), DIMENSION(:) :: qp_eigval
664 : INTEGER :: ispin
665 :
666 : CHARACTER(LEN=*), PARAMETER :: routineN = 'print_static_cohsex_mos'
667 :
668 : CHARACTER(LEN=5) :: spin_label
669 : CHARACTER(LEN=default_path_length) :: fname, title
670 : INTEGER :: first_mo, handle, ivector, last_mo, &
671 : n_mo, nelectron, nhomo, nlumo, unit_nr
672 6 : INTEGER, DIMENSION(:), POINTER :: stride
673 : LOGICAL :: print_coeff, print_cubes, &
674 : print_energies, print_mo, print_occ
675 : REAL(KIND=dp) :: maxocc
676 6 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
677 : TYPE(cell_type), POINTER :: cell
678 : TYPE(cp_logger_type), POINTER :: logger
679 : TYPE(dft_control_type), POINTER :: dft_control
680 : TYPE(mo_set_type) :: qp_mo_set
681 : TYPE(particle_list_type), POINTER :: particles
682 6 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
683 : TYPE(pw_c1d_gs_type) :: wf_g
684 : TYPE(pw_env_type), POINTER :: pw_env
685 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
686 : TYPE(pw_r3d_rs_type) :: wf_r
687 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
688 : TYPE(qs_subsys_type), POINTER :: subsys
689 : TYPE(section_vals_type), POINTER :: dft_section, input
690 :
691 6 : CALL timeset(routineN, handle)
692 :
693 6 : n_mo = bs_env%n_ao
694 :
695 6 : CALL get_qs_env(qs_env, input=input)
696 6 : dft_section => section_vals_get_subs_vals(input, "DFT")
697 6 : logger => cp_get_default_logger()
698 :
699 : ! ---------------------------------------------------------------------------------------
700 : ! MO energies / occupation numbers / coefficients
701 : ! ---------------------------------------------------------------------------------------
702 6 : CALL section_vals_val_get(dft_section, "PRINT%MO%ENERGIES", l_val=print_energies)
703 6 : CALL section_vals_val_get(dft_section, "PRINT%MO%COEFFICIENTS", l_val=print_coeff)
704 6 : CALL section_vals_val_get(dft_section, "PRINT%MO%OCCUPATION_NUMBERS", l_val=print_occ)
705 6 : print_mo = (cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO") /= 0)
706 :
707 6 : IF (print_mo .AND. (print_energies .OR. print_occ .OR. print_coeff)) THEN
708 :
709 : ! occupation numbers of the QP states: the lowest n_occ are occupied
710 0 : IF (bs_env%n_spin == 1) THEN
711 0 : maxocc = 2.0_dp
712 : ELSE
713 0 : maxocc = 1.0_dp
714 : END IF
715 0 : nelectron = NINT(maxocc*REAL(bs_env%n_occ(ispin), KIND=dp))
716 :
717 : CALL allocate_mo_set(qp_mo_set, nao=n_mo, nmo=n_mo, nelectron=nelectron, &
718 : n_el_f=REAL(nelectron, KIND=dp), maxocc=maxocc, &
719 0 : flexible_electron_count=0.0_dp)
720 0 : CALL init_mo_set(qp_mo_set, fm_struct=fm_qp_mo_ao%matrix_struct, name="static_cohsex")
721 0 : qp_mo_set%eigenvalues(:) = qp_eigval(:)
722 0 : CALL cp_fm_to_fm(fm_qp_mo_ao, qp_mo_set%mo_coeff)
723 :
724 0 : CALL get_qs_env(qs_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
725 :
726 0 : IF (bs_env%n_spin == 1) THEN
727 : CALL write_mo_set_to_output_unit(qp_mo_set, qs_kind_set, particle_set, dft_section, &
728 0 : 4, 0, final_mos=.TRUE.)
729 : ELSE
730 0 : IF (ispin == 1) THEN
731 0 : spin_label = "ALPHA"
732 : ELSE
733 0 : spin_label = "BETA"
734 : END IF
735 : CALL write_mo_set_to_output_unit(qp_mo_set, qs_kind_set, particle_set, dft_section, &
736 0 : 4, 0, final_mos=.TRUE., spin=TRIM(spin_label))
737 : END IF
738 :
739 0 : CALL deallocate_mo_set(qp_mo_set)
740 : END IF
741 :
742 : ! ---------------------------------------------------------------------------------------
743 : ! MO cubes
744 : ! ---------------------------------------------------------------------------------------
745 6 : print_cubes = (cp_print_key_should_output(logger%iter_info, dft_section, "PRINT%MO_CUBES") /= 0)
746 :
747 6 : IF (print_cubes) THEN
748 0 : CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%NHOMO", i_val=nhomo)
749 0 : CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%NLUMO", i_val=nlumo)
750 0 : CALL section_vals_val_get(dft_section, "PRINT%MO_CUBES%STRIDE", i_vals=stride)
751 :
752 0 : IF (nhomo == -1) nhomo = bs_env%n_occ(ispin)
753 0 : IF (nlumo == -1) nlumo = n_mo - bs_env%n_occ(ispin)
754 0 : first_mo = MAX(1, bs_env%n_occ(ispin) - nhomo + 1)
755 0 : last_mo = MIN(n_mo, bs_env%n_occ(ispin) + nlumo)
756 :
757 : CALL get_qs_env(qs_env, &
758 : atomic_kind_set=atomic_kind_set, &
759 : qs_kind_set=qs_kind_set, &
760 : cell=cell, &
761 : dft_control=dft_control, &
762 : particle_set=particle_set, &
763 : pw_env=pw_env, &
764 0 : subsys=subsys)
765 0 : CALL qs_subsys_get(subsys, particles=particles)
766 :
767 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
768 0 : CALL auxbas_pw_pool%create_pw(wf_r)
769 0 : CALL auxbas_pw_pool%create_pw(wf_g)
770 :
771 0 : DO ivector = first_mo, last_mo
772 : CALL calculate_wavefunction(fm_qp_mo_ao, ivector, wf_r, wf_g, atomic_kind_set, &
773 0 : qs_kind_set, cell, dft_control, particle_set, pw_env)
774 :
775 0 : WRITE (fname, '(A,I5.5,A,I1)') "static_cohsex_WFN_", ivector, "_", ispin
776 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%MO_CUBES", &
777 : extension=".cube", middle_name=TRIM(fname), &
778 0 : file_position="REWIND", log_filename=.FALSE.)
779 0 : WRITE (title, *) "static COHSEX WAVEFUNCTION ", ivector, " spin ", ispin, &
780 0 : " i.e. HOMO ", ivector - bs_env%n_occ(ispin)
781 0 : CALL cp_pw_to_cube(wf_r, unit_nr, title, particles, stride=stride)
782 0 : CALL cp_print_key_finished_output(unit_nr, logger, dft_section, "PRINT%MO_CUBES")
783 : END DO
784 :
785 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
786 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
787 : END IF
788 :
789 6 : CALL timestop(handle)
790 :
791 6 : END SUBROUTINE print_static_cohsex_mos
792 :
793 : ! **************************************************************************************************
794 : !> \brief ...
795 : !> \param fm ...
796 : !> \param matrix_index ...
797 : !> \param matrix_name ...
798 : !> \param qs_env ...
799 : ! **************************************************************************************************
800 6 : SUBROUTINE fm_write_formatted(fm, matrix_index, matrix_name, qs_env)
801 : TYPE(cp_fm_type) :: fm
802 : INTEGER :: matrix_index
803 : CHARACTER(LEN=*) :: matrix_name
804 : TYPE(qs_environment_type), POINTER :: qs_env
805 :
806 : CHARACTER(LEN=*), PARAMETER :: key = 'PROPERTIES%BANDSTRUCTURE%GW%PRINT%RESTART', &
807 : routineN = 'fm_write_formatted'
808 :
809 : CHARACTER(LEN=default_string_length) :: filename
810 : INTEGER :: handle, unit_nr
811 : TYPE(cp_logger_type), POINTER :: logger
812 : TYPE(section_vals_type), POINTER :: input
813 :
814 6 : CALL timeset(routineN, handle)
815 :
816 6 : CALL get_qs_env(qs_env, input=input)
817 :
818 6 : logger => cp_get_default_logger()
819 :
820 6 : IF (BTEST(cp_print_key_should_output(logger%iter_info, input, key), cp_p_file)) THEN
821 :
822 0 : IF (matrix_index < 10) THEN
823 0 : WRITE (filename, '(3A,I1)') "RESTART_", matrix_name, "_0", matrix_index
824 0 : ELSE IF (matrix_index < 100) THEN
825 0 : WRITE (filename, '(3A,I2)') "RESTART_", matrix_name, "_", matrix_index
826 : ELSE
827 0 : CPABORT('Please implement more than 99 time/frequency points.')
828 : END IF
829 :
830 : unit_nr = cp_print_key_unit_nr(logger, input, key, extension=".matrix", &
831 : file_form="FORMATTED", middle_name=TRIM(filename), &
832 0 : file_position="REWIND", file_action="WRITE")
833 :
834 0 : CALL cp_fm_write_formatted(fm, unit_nr)
835 0 : IF (unit_nr > 0) THEN
836 0 : CALL close_file(unit_nr)
837 : END IF
838 : END IF
839 :
840 6 : CALL timestop(handle)
841 :
842 6 : END SUBROUTINE fm_write_formatted
843 :
844 : END MODULE gw_tensor_large_cell_Gamma_static_cohsex
|
