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 Does all kind of post scf calculations for GPW/GAPW
10 : !> \par History
11 : !> Started as a copy from the relevant part of qs_scf
12 : !> Start to adapt for k-points [07.2015, JGH]
13 : !> \author Joost VandeVondele (10.2003)
14 : ! **************************************************************************************************
15 : MODULE qs_energy_window
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_dbcsr_api, ONLY: &
18 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_get_info, dbcsr_multiply, &
19 : dbcsr_p_type, dbcsr_release, dbcsr_type
20 : USE cp_dbcsr_contrib, ONLY: dbcsr_frobenius_norm
21 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
22 : copy_fm_to_dbcsr
23 : USE cp_fm_basic_linalg, ONLY: cp_fm_trace
24 : USE cp_fm_diag, ONLY: choose_eigv_solver
25 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
26 : cp_fm_struct_release,&
27 : cp_fm_struct_type
28 : USE cp_fm_types, ONLY: cp_fm_create,&
29 : cp_fm_release,&
30 : cp_fm_to_fm,&
31 : cp_fm_type
32 : USE cp_log_handling, ONLY: cp_get_default_logger,&
33 : cp_logger_get_default_io_unit,&
34 : cp_logger_type
35 : USE cp_output_handling, ONLY: cp_iter_string,&
36 : cp_print_key_finished_output,&
37 : cp_print_key_unit_nr
38 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
39 : USE input_section_types, ONLY: section_get_ivals,&
40 : section_vals_get_subs_vals,&
41 : section_vals_type,&
42 : section_vals_val_get
43 : USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz
44 : USE kinds, ONLY: dp
45 : USE parallel_gemm_api, ONLY: parallel_gemm
46 : USE particle_list_types, ONLY: particle_list_type
47 : USE pw_env_types, ONLY: pw_env_get,&
48 : pw_env_type
49 : USE pw_methods, ONLY: pw_integrate_function
50 : USE pw_pool_types, ONLY: pw_pool_type
51 : USE pw_types, ONLY: pw_c1d_gs_type,&
52 : pw_r3d_rs_type
53 : USE qs_collocate_density, ONLY: calculate_rho_elec
54 : USE qs_environment_types, ONLY: get_qs_env,&
55 : qs_environment_type
56 : USE qs_ks_types, ONLY: qs_ks_env_type
57 : USE qs_rho_types, ONLY: qs_rho_get,&
58 : qs_rho_type
59 : USE qs_subsys_types, ONLY: qs_subsys_get,&
60 : qs_subsys_type
61 : #include "./base/base_uses.f90"
62 :
63 : IMPLICIT NONE
64 : PRIVATE
65 :
66 : ! Global parameters
67 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_energy_window'
68 :
69 : PUBLIC :: energy_windows
70 :
71 : ! **************************************************************************************************
72 :
73 : CONTAINS
74 :
75 : ! **************************************************************************************************
76 : !> \brief ...
77 : !> \param qs_env ...
78 : ! **************************************************************************************************
79 112 : SUBROUTINE energy_windows(qs_env)
80 :
81 : TYPE(qs_environment_type), POINTER :: qs_env
82 :
83 : CHARACTER(len=*), PARAMETER :: routineN = 'energy_windows'
84 : LOGICAL, PARAMETER :: debug = .FALSE.
85 : REAL(KIND=dp), PARAMETER :: one = 1.0_dp, zero = 0.0_dp
86 :
87 : CHARACTER(len=40) :: ext, title
88 : INTEGER :: handle, i, lanzcos_max_iter, last, nao, &
89 : nelectron_total, newton_schulz_order, &
90 : next, nwindows, print_unit, unit_nr
91 112 : INTEGER, DIMENSION(:), POINTER :: stride(:)
92 : LOGICAL :: mpi_io, print_cube, restrict_range
93 : REAL(KIND=dp) :: bin_width, density_ewindow_total, density_total, energy_range, fermi_level, &
94 : filter_eps, frob_norm, lanzcos_threshold, lower_bound, occupation, upper_bound
95 112 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, P_eigenvalues, &
96 112 : window_eigenvalues
97 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
98 : TYPE(cp_fm_struct_type), POINTER :: ao_ao_fmstruct, window_fm_struct
99 : TYPE(cp_fm_type) :: eigenvectors, eigenvectors_nonorth, matrix_ks_fm, P_eigenvectors, &
100 : P_window_fm, rho_ao_ortho_fm, S_minus_half_fm, tmp_fm, window_eigenvectors, window_fm
101 : TYPE(cp_logger_type), POINTER :: logger
102 112 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s, rho_ao
103 : TYPE(dbcsr_type) :: matrix_ks_nosym, S_half, S_minus_half, &
104 : tmp
105 : TYPE(dbcsr_type), POINTER :: rho_ao_ortho, window
106 : TYPE(particle_list_type), POINTER :: particles
107 : TYPE(pw_c1d_gs_type) :: rho_g
108 : TYPE(pw_env_type), POINTER :: pw_env
109 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
110 : TYPE(pw_r3d_rs_type) :: rho_r
111 : TYPE(qs_ks_env_type), POINTER :: ks_env
112 : TYPE(qs_rho_type), POINTER :: rho
113 : TYPE(qs_subsys_type), POINTER :: subsys
114 : TYPE(section_vals_type), POINTER :: dft_section, input, ls_scf_section
115 :
116 112 : CALL timeset(routineN, handle)
117 :
118 112 : logger => cp_get_default_logger()
119 112 : unit_nr = cp_logger_get_default_io_unit(logger)
120 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, matrix_ks=matrix_ks, pw_env=pw_env, rho=rho, &
121 112 : input=input, nelectron_total=nelectron_total, subsys=subsys, ks_env=ks_env, matrix_s=matrix_s)
122 112 : CALL qs_subsys_get(subsys, particles=particles)
123 112 : CALL qs_rho_get(rho_struct=rho, rho_ao=rho_ao)
124 112 : IF (SIZE(rho_ao) > 1) CALL cp_warn(__LOCATION__, &
125 0 : "The printing of energy windows is only implemented for closed shell systems")
126 112 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
127 :
128 : !reading the input
129 112 : dft_section => section_vals_get_subs_vals(input, "DFT")
130 112 : ls_scf_section => section_vals_get_subs_vals(input, "DFT%LS_SCF")
131 112 : CALL section_vals_val_get(dft_section, "PRINT%ENERGY_WINDOWS%N_WINDOWS", i_val=nwindows)
132 112 : CALL section_vals_val_get(dft_section, "PRINT%ENERGY_WINDOWS%PRINT_CUBES", l_val=print_cube)
133 112 : CALL section_vals_val_get(dft_section, "PRINT%ENERGY_WINDOWS%RESTRICT_RANGE", l_val=restrict_range)
134 112 : CALL section_vals_val_get(dft_section, "PRINT%ENERGY_WINDOWS%RANGE", r_val=energy_range)
135 : NULLIFY (stride)
136 112 : ALLOCATE (stride(3))
137 448 : stride = section_get_ivals(dft_section, "PRINT%ENERGY_WINDOWS%STRIDE")
138 112 : CALL section_vals_val_get(dft_section, "PRINT%ENERGY_WINDOWS%EPS_FILTER", r_val=filter_eps)
139 112 : CALL section_vals_val_get(ls_scf_section, "EPS_LANCZOS", r_val=lanzcos_threshold)
140 112 : CALL section_vals_val_get(ls_scf_section, "MAX_ITER_LANCZOS", i_val=lanzcos_max_iter)
141 112 : CALL section_vals_val_get(ls_scf_section, "SIGN_SQRT_ORDER", i_val=newton_schulz_order)
142 :
143 : !Initialize data
144 112 : ALLOCATE (window, rho_ao_ortho)
145 112 : CALL dbcsr_get_info(matrix=matrix_ks(1)%matrix, nfullrows_total=nao)
146 336 : ALLOCATE (eigenvalues(nao))
147 112 : CALL dbcsr_create(tmp, template=matrix_s(1)%matrix, matrix_type="N")
148 112 : CALL dbcsr_create(S_minus_half, template=matrix_s(1)%matrix, matrix_type="N")
149 112 : CALL dbcsr_create(S_half, template=matrix_s(1)%matrix, matrix_type="N")
150 112 : CALL dbcsr_create(window, template=matrix_s(1)%matrix, matrix_type="N")
151 112 : CALL dbcsr_create(rho_ao_ortho, template=matrix_s(1)%matrix, matrix_type="N")
152 112 : CALL cp_fm_struct_create(fmstruct=ao_ao_fmstruct, context=blacs_env, nrow_global=nao, ncol_global=nao)
153 112 : CALL cp_fm_create(P_window_fm, ao_ao_fmstruct)
154 112 : CALL cp_fm_create(matrix_ks_fm, ao_ao_fmstruct)
155 112 : CALL cp_fm_create(rho_ao_ortho_fm, ao_ao_fmstruct)
156 112 : CALL cp_fm_create(S_minus_half_fm, ao_ao_fmstruct)
157 112 : CALL cp_fm_create(eigenvectors, ao_ao_fmstruct)
158 112 : CALL cp_fm_create(eigenvectors_nonorth, ao_ao_fmstruct)
159 112 : CALL auxbas_pw_pool%create_pw(rho_r)
160 112 : CALL auxbas_pw_pool%create_pw(rho_g)
161 :
162 : !calculate S_minus_half
163 : CALL matrix_sqrt_Newton_Schulz(S_half, S_minus_half, matrix_s(1)%matrix, filter_eps, &
164 112 : newton_schulz_order, lanzcos_threshold, lanzcos_max_iter)
165 :
166 : !get the full ks matrix
167 112 : CALL dbcsr_desymmetrize(matrix_ks(1)%matrix, matrix_ks_nosym)
168 :
169 : !switching to orthonormal basis
170 112 : CALL dbcsr_multiply("N", "N", one, S_minus_half, matrix_ks_nosym, zero, tmp, filter_eps=filter_eps)
171 112 : CALL dbcsr_multiply("N", "N", one, tmp, S_minus_half, zero, matrix_ks_nosym, filter_eps=filter_eps)
172 112 : CALL copy_dbcsr_to_fm(matrix_ks_nosym, matrix_ks_fm)
173 112 : CALL dbcsr_multiply("N", "N", one, S_half, rho_ao(1)%matrix, zero, tmp, filter_eps=filter_eps)
174 112 : CALL dbcsr_multiply("N", "N", one, tmp, S_half, zero, rho_ao_ortho, filter_eps=filter_eps)
175 112 : CALL copy_dbcsr_to_fm(rho_ao_ortho, rho_ao_ortho_fm)
176 :
177 : !diagonalize the full ks matrix
178 112 : CALL choose_eigv_solver(matrix_ks_fm, eigenvectors, eigenvalues)
179 112 : fermi_level = eigenvalues((nelectron_total + MOD(nelectron_total, 2))/2)
180 112 : IF (restrict_range) THEN
181 12 : lower_bound = MAX(fermi_level - energy_range, eigenvalues(1))
182 12 : upper_bound = MIN(fermi_level + energy_range, eigenvalues(SIZE(eigenvalues)))
183 : ELSE
184 100 : lower_bound = eigenvalues(1)
185 100 : upper_bound = eigenvalues(SIZE(eigenvalues))
186 : END IF
187 112 : IF (unit_nr > 0) THEN
188 56 : WRITE (unit_nr, *) " Creating energy windows. Fermi level: ", fermi_level
189 56 : WRITE (unit_nr, *) " Printing Energy Levels from ", lower_bound, " to ", upper_bound
190 : END IF
191 : !Rotate the eigenvectors back out of the orthonormal basis
192 112 : CALL copy_dbcsr_to_fm(S_minus_half, S_minus_half_fm)
193 : !calculate the density caused by the mos in the energy window
194 112 : CALL parallel_gemm("N", "N", nao, nao, nao, one, S_minus_half_fm, eigenvectors, zero, eigenvectors_nonorth)
195 :
196 : IF (debug) THEN
197 : !check difference to actual density
198 : CALL cp_fm_struct_create(fmstruct=window_fm_struct, context=blacs_env, nrow_global=nao, &
199 : ncol_global=nelectron_total/2)
200 : CALL cp_fm_create(window_fm, window_fm_struct)
201 : CALL cp_fm_to_fm(eigenvectors_nonorth, window_fm, nelectron_total/2, 1, 1)
202 : CALL parallel_gemm("N", "T", nao, nao, nelectron_total/2, 2*one, window_fm, window_fm, zero, P_window_fm)
203 : !ensure the correct sparsity
204 : CALL copy_fm_to_dbcsr(P_window_fm, tmp)
205 : CALL dbcsr_copy(window, matrix_ks(1)%matrix)
206 : CALL dbcsr_copy(window, tmp, keep_sparsity=.TRUE.)
207 : CALL calculate_rho_elec(matrix_p=window, &
208 : rho=rho_r, &
209 : rho_gspace=rho_g, &
210 : ks_env=ks_env)
211 : density_total = pw_integrate_function(rho_r)
212 : IF (unit_nr > 0) WRITE (unit_nr, *) " Ground-state density: ", density_total
213 : frob_norm = dbcsr_frobenius_norm(window)
214 : IF (unit_nr > 0) WRITE (unit_nr, *) " Frob norm of calculated ground-state density matrix: ", frob_norm
215 : CALL dbcsr_add(window, rho_ao(1)%matrix, one, -one)
216 : frob_norm = dbcsr_frobenius_norm(rho_ao(1)%matrix)
217 : IF (unit_nr > 0) WRITE (unit_nr, *) " Frob norm of current density matrix: ", frob_norm
218 : frob_norm = dbcsr_frobenius_norm(window)
219 : IF (unit_nr > 0) WRITE (unit_nr, *) " Difference between calculated ground-state density and current density: ", frob_norm
220 : CALL cp_fm_struct_release(window_fm_struct)
221 : CALL cp_fm_create(tmp_fm, ao_ao_fmstruct)
222 : CALL cp_fm_to_fm(rho_ao_ortho_fm, tmp_fm)
223 : CALL cp_fm_create(P_eigenvectors, ao_ao_fmstruct)
224 : ALLOCATE (P_eigenvalues(nao))
225 : CALL choose_eigv_solver(tmp_fm, P_eigenvectors, P_eigenvalues)
226 : CALL cp_fm_create(window_eigenvectors, ao_ao_fmstruct)
227 : ALLOCATE (window_eigenvalues(nao))
228 : CALL cp_fm_to_fm(eigenvectors, window_fm, nelectron_total/2, 1, 1)
229 : CALL parallel_gemm("N", "T", nao, nao, nelectron_total/2, 2*one, window_fm, window_fm, zero, P_window_fm)
230 : CALL choose_eigv_solver(P_window_fm, window_eigenvectors, window_eigenvalues)
231 : DO i = 1, nao
232 : IF (unit_nr > 0) THEN
233 : WRITE (unit_nr, *) i, "H:", eigenvalues(i), "P:", P_eigenvalues(nao - i + 1), "Pnew:", window_eigenvalues(nao - i + 1)
234 : END IF
235 : END DO
236 : DEALLOCATE (P_eigenvalues)
237 : CALL cp_fm_release(tmp_fm)
238 : CALL cp_fm_release(P_eigenvectors)
239 : DEALLOCATE (window_eigenvalues)
240 : CALL cp_fm_release(window_eigenvectors)
241 : CALL cp_fm_release(window_fm)
242 : END IF
243 :
244 : !create energy windows
245 112 : bin_width = (upper_bound - lower_bound)/nwindows
246 112 : next = 0
247 :
248 9312 : DO i = 1, nwindows
249 9210 : DO WHILE (eigenvalues(next + 1) < lower_bound)
250 9200 : next = next + 1
251 : END DO
252 : last = next
253 11724 : DO WHILE (eigenvalues(next + 1) < lower_bound + i*bin_width)
254 2544 : next = next + 1
255 11724 : IF (next == SIZE(eigenvalues)) EXIT
256 : END DO
257 : !calculate the occupation
258 : !not sure how bad this is now load balanced due to using the same blacs_env
259 9200 : CALL cp_fm_struct_create(fmstruct=window_fm_struct, context=blacs_env, nrow_global=nao, ncol_global=next - last)
260 9200 : CALL cp_fm_create(window_fm, window_fm_struct)
261 : !copy the mos in the energy window into a separate matrix
262 9200 : CALL cp_fm_to_fm(eigenvectors, window_fm, next - last, last + 1, 1)
263 9200 : CALL parallel_gemm("N", "T", nao, nao, next - last, one, window_fm, window_fm, zero, P_window_fm)
264 9200 : CALL cp_fm_trace(P_window_fm, rho_ao_ortho_fm, occupation)
265 9200 : IF (print_cube) THEN
266 180 : CALL cp_fm_to_fm(eigenvectors_nonorth, window_fm, next - last, last + 1, 1)
267 : !print the energy window to a cube file
268 : !calculate the density caused by the mos in the energy window
269 180 : CALL parallel_gemm("N", "T", nao, nao, next - last, one, window_fm, window_fm, zero, P_window_fm)
270 180 : CALL copy_fm_to_dbcsr(P_window_fm, tmp)
271 : !ensure the correct sparsity
272 180 : CALL dbcsr_copy(window, matrix_ks(1)%matrix)
273 180 : CALL dbcsr_copy(window, tmp, keep_sparsity=.TRUE.)
274 : CALL calculate_rho_elec(matrix_p=window, &
275 : rho=rho_r, &
276 : rho_gspace=rho_g, &
277 180 : ks_env=ks_env)
278 180 : WRITE (ext, "(A14,I5.5,A)") "-ENERGY-WINDOW", i, TRIM(cp_iter_string(logger%iter_info))//".cube"
279 180 : mpi_io = .TRUE.
280 : print_unit = cp_print_key_unit_nr(logger, dft_section, "PRINT%ENERGY_WINDOWS", &
281 : extension=ext, file_status="REPLACE", file_action="WRITE", &
282 180 : log_filename=.FALSE., mpi_io=mpi_io)
283 180 : WRITE (title, "(A14,I5)") "ENERGY WINDOW ", i
284 180 : CALL cp_pw_to_cube(rho_r, print_unit, title, particles=particles, stride=stride, mpi_io=mpi_io)
285 : CALL cp_print_key_finished_output(print_unit, logger, dft_section, &
286 180 : "PRINT%ENERGY_WINDOWS", mpi_io=mpi_io)
287 180 : density_ewindow_total = pw_integrate_function(rho_r)
288 270 : IF (unit_nr > 0) WRITE (unit_nr, "(A,F16.10,A,I5,A,F20.14,A,F20.14)") " Energy Level: ", &
289 90 : lower_bound + (i - 0.5_dp)*bin_width, " Number of states: ", next - last, " Occupation: ", &
290 180 : occupation, " Grid Density ", density_ewindow_total
291 : ELSE
292 9020 : IF (unit_nr > 0) THEN
293 4510 : WRITE (unit_nr, "(A,F16.10,A,I5,A,F20.14)") " Energy Level: ", lower_bound + (i - 0.5_dp)*bin_width, &
294 9020 : " Number of states: ", next - last, " Occupation: ", occupation
295 : END IF
296 : END IF
297 9200 : CALL cp_fm_release(window_fm)
298 27712 : CALL cp_fm_struct_release(window_fm_struct)
299 : END DO
300 :
301 : !clean up
302 112 : CALL dbcsr_release(matrix_ks_nosym)
303 112 : CALL dbcsr_release(tmp)
304 112 : CALL dbcsr_release(window)
305 112 : CALL dbcsr_release(S_minus_half)
306 112 : CALL dbcsr_release(S_half)
307 112 : CALL dbcsr_release(rho_ao_ortho)
308 112 : DEALLOCATE (window, rho_ao_ortho)
309 112 : CALL cp_fm_struct_release(ao_ao_fmstruct)
310 112 : CALL cp_fm_release(matrix_ks_fm)
311 112 : CALL cp_fm_release(rho_ao_ortho_fm)
312 112 : CALL cp_fm_release(eigenvectors)
313 112 : CALL cp_fm_release(P_window_fm)
314 112 : CALL cp_fm_release(eigenvectors_nonorth)
315 112 : CALL cp_fm_release(S_minus_half_fm)
316 112 : CALL auxbas_pw_pool%give_back_pw(rho_r)
317 112 : CALL auxbas_pw_pool%give_back_pw(rho_g)
318 112 : DEALLOCATE (eigenvalues)
319 112 : DEALLOCATE (STRIDE)
320 :
321 112 : CALL timestop(handle)
322 :
323 560 : END SUBROUTINE energy_windows
324 :
325 : !**************************************************************************************************
326 :
327 : END MODULE qs_energy_window
|
