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