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