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 Routines for a linear scaling quickstep SCF run based on the density
10 : !> matrix
11 : !> \par History
12 : !> 2010.10 created [Joost VandeVondele]
13 : !> \author Joost VandeVondele
14 : ! **************************************************************************************************
15 : MODULE dm_ls_scf
16 : USE arnoldi_api, ONLY: arnoldi_extremal
17 : USE bibliography, ONLY: Kolafa2004,&
18 : Kuhne2007,&
19 : cite_reference
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_dbcsr_api, ONLY: &
22 : dbcsr_add, dbcsr_binary_read, dbcsr_binary_write, dbcsr_copy, dbcsr_create, &
23 : dbcsr_distribution_type, dbcsr_filter, dbcsr_get_info, dbcsr_get_occupation, &
24 : dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_type, &
25 : dbcsr_type_no_symmetry
26 : USE cp_dbcsr_contrib, ONLY: dbcsr_checksum
27 : USE cp_external_control, ONLY: external_control
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_get_default_unit_nr,&
30 : cp_logger_type
31 : USE dm_ls_chebyshev, ONLY: compute_chebyshev
32 : USE dm_ls_scf_create, ONLY: ls_scf_create
33 : USE dm_ls_scf_curvy, ONLY: deallocate_curvy_data,&
34 : dm_ls_curvy_optimization
35 : USE dm_ls_scf_methods, ONLY: apply_matrix_preconditioner,&
36 : compute_homo_lumo,&
37 : density_matrix_sign,&
38 : density_matrix_sign_fixed_mu,&
39 : density_matrix_tc2,&
40 : density_matrix_trs4,&
41 : ls_scf_init_matrix_S
42 : USE dm_ls_scf_qs, ONLY: &
43 : ls_nonscf_energy, ls_nonscf_ks, ls_scf_dm_to_ks, ls_scf_init_qs, ls_scf_qs_atomic_guess, &
44 : matrix_ls_create, matrix_ls_to_qs, matrix_qs_to_ls, rho_mixing_ls_init
45 : USE dm_ls_scf_types, ONLY: ls_scf_env_type
46 : USE ec_env_types, ONLY: energy_correction_type
47 : USE input_constants, ONLY: ls_cluster_atomic,&
48 : ls_scf_pexsi,&
49 : ls_scf_sign,&
50 : ls_scf_tc2,&
51 : ls_scf_trs4,&
52 : transport_transmission
53 : USE input_section_types, ONLY: section_vals_type
54 : USE iterate_matrix, ONLY: purify_mcweeny
55 : USE kinds, ONLY: default_path_length,&
56 : default_string_length,&
57 : dp
58 : USE machine, ONLY: m_flush,&
59 : m_walltime
60 : USE mathlib, ONLY: binomial
61 : USE molecule_types, ONLY: molecule_type
62 : USE pao_main, ONLY: pao_optimization_end,&
63 : pao_optimization_start,&
64 : pao_post_scf,&
65 : pao_update
66 : USE pexsi_methods, ONLY: density_matrix_pexsi,&
67 : pexsi_finalize_scf,&
68 : pexsi_init_scf,&
69 : pexsi_set_convergence_tolerance,&
70 : pexsi_to_qs
71 : USE qs_diis, ONLY: qs_diis_b_clear_sparse,&
72 : qs_diis_b_create_sparse,&
73 : qs_diis_b_step_4lscf
74 : USE qs_diis_types, ONLY: qs_diis_b_release_sparse,&
75 : qs_diis_buffer_type_sparse
76 : USE qs_environment_types, ONLY: get_qs_env,&
77 : qs_environment_type
78 : USE qs_ks_types, ONLY: qs_ks_env_type
79 : USE qs_nonscf_utils, ONLY: qs_nonscf_print_summary
80 : USE qs_scf_post_gpw, ONLY: qs_scf_post_moments,&
81 : write_mo_free_results
82 : USE qs_scf_post_tb, ONLY: scf_post_calculation_tb
83 : USE transport, ONLY: external_scf_method,&
84 : transport_initialize
85 : USE transport_env_types, ONLY: transport_env_type
86 : #include "./base/base_uses.f90"
87 :
88 : IMPLICIT NONE
89 :
90 : PRIVATE
91 :
92 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf'
93 :
94 : PUBLIC :: calculate_w_matrix_ls, ls_scf, post_scf_sparsities
95 :
96 : CONTAINS
97 :
98 : ! **************************************************************************************************
99 : !> \brief perform an linear scaling scf procedure: entry point
100 : !>
101 : !> \param qs_env ...
102 : !> \param nonscf ...
103 : !> \par History
104 : !> 2010.10 created [Joost VandeVondele]
105 : !> \author Joost VandeVondele
106 : ! **************************************************************************************************
107 664 : SUBROUTINE ls_scf(qs_env, nonscf)
108 : TYPE(qs_environment_type), POINTER :: qs_env
109 : LOGICAL, INTENT(IN), OPTIONAL :: nonscf
110 :
111 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf'
112 :
113 : INTEGER :: handle
114 : LOGICAL :: do_scf, pao_is_done
115 : TYPE(ls_scf_env_type), POINTER :: ls_scf_env
116 :
117 664 : CALL timeset(routineN, handle)
118 664 : do_scf = .TRUE.
119 664 : IF (PRESENT(nonscf)) do_scf = .NOT. nonscf
120 :
121 : ! Moved here from qs_environment to remove dependencies
122 664 : CALL ls_scf_create(qs_env)
123 664 : CALL get_qs_env(qs_env, ls_scf_env=ls_scf_env)
124 :
125 664 : IF (do_scf) THEN
126 598 : CALL pao_optimization_start(qs_env, ls_scf_env)
127 598 : pao_is_done = .FALSE.
128 1414 : DO WHILE (.NOT. pao_is_done)
129 816 : CALL ls_scf_init_scf(qs_env, ls_scf_env, .FALSE.)
130 816 : CALL pao_update(qs_env, ls_scf_env, pao_is_done)
131 816 : CALL ls_scf_main(qs_env, ls_scf_env, .FALSE.)
132 816 : CALL pao_post_scf(qs_env, ls_scf_env, pao_is_done)
133 816 : CALL ls_scf_post(qs_env, ls_scf_env)
134 : END DO
135 598 : CALL pao_optimization_end(ls_scf_env)
136 : ELSE
137 66 : CALL ls_scf_init_scf(qs_env, ls_scf_env, .TRUE.)
138 66 : CALL ls_scf_main(qs_env, ls_scf_env, .TRUE.)
139 66 : CALL ls_scf_post(qs_env, ls_scf_env)
140 : END IF
141 :
142 664 : CALL timestop(handle)
143 :
144 664 : END SUBROUTINE ls_scf
145 :
146 : ! **************************************************************************************************
147 : !> \brief initialization needed for scf
148 : !> \param qs_env ...
149 : !> \param ls_scf_env ...
150 : !> \param nonscf ...
151 : !> \par History
152 : !> 2010.10 created [Joost VandeVondele]
153 : !> \author Joost VandeVondele
154 : ! **************************************************************************************************
155 882 : SUBROUTINE ls_scf_init_scf(qs_env, ls_scf_env, nonscf)
156 : TYPE(qs_environment_type), POINTER :: qs_env
157 : TYPE(ls_scf_env_type) :: ls_scf_env
158 : LOGICAL, INTENT(IN) :: nonscf
159 :
160 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_init_scf'
161 :
162 : INTEGER :: handle, ispin, nspin, unit_nr
163 : TYPE(cp_logger_type), POINTER :: logger
164 882 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_w
165 : TYPE(dft_control_type), POINTER :: dft_control
166 882 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
167 : TYPE(qs_ks_env_type), POINTER :: ks_env
168 : TYPE(section_vals_type), POINTER :: input
169 :
170 882 : CALL timeset(routineN, handle)
171 :
172 : ! get a useful output_unit
173 882 : logger => cp_get_default_logger()
174 882 : IF (logger%para_env%is_source()) THEN
175 441 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
176 : ELSE
177 : unit_nr = -1
178 : END IF
179 :
180 : ! get basic quantities from the qs_env
181 : CALL get_qs_env(qs_env, nelectron_total=ls_scf_env%nelectron_total, &
182 : matrix_s=matrix_s, &
183 : matrix_w=matrix_w, &
184 : ks_env=ks_env, &
185 : dft_control=dft_control, &
186 : molecule_set=molecule_set, &
187 : input=input, &
188 : has_unit_metric=ls_scf_env%has_unit_metric, &
189 : para_env=ls_scf_env%para_env, &
190 882 : nelectron_spin=ls_scf_env%nelectron_spin)
191 :
192 : ! needs forces ? There might be a better way to flag this
193 882 : ls_scf_env%calculate_forces = ASSOCIATED(matrix_w)
194 :
195 : ! some basic initialization of the QS side of things
196 882 : CALL ls_scf_init_qs(qs_env)
197 :
198 : ! create the matrix template for use in the ls procedures
199 : CALL matrix_ls_create(matrix_ls=ls_scf_env%matrix_s, matrix_qs=matrix_s(1)%matrix, &
200 882 : ls_mstruct=ls_scf_env%ls_mstruct)
201 :
202 882 : nspin = ls_scf_env%nspins
203 882 : IF (ALLOCATED(ls_scf_env%matrix_p)) THEN
204 1098 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
205 1098 : CALL dbcsr_release(ls_scf_env%matrix_p(ispin))
206 : END DO
207 : ELSE
208 1358 : ALLOCATE (ls_scf_env%matrix_p(nspin))
209 : END IF
210 :
211 1784 : DO ispin = 1, nspin
212 : CALL dbcsr_create(ls_scf_env%matrix_p(ispin), template=ls_scf_env%matrix_s, &
213 1784 : matrix_type=dbcsr_type_no_symmetry)
214 : END DO
215 :
216 3548 : ALLOCATE (ls_scf_env%matrix_ks(nspin))
217 1784 : DO ispin = 1, nspin
218 : CALL dbcsr_create(ls_scf_env%matrix_ks(ispin), template=ls_scf_env%matrix_s, &
219 1784 : matrix_type=dbcsr_type_no_symmetry)
220 : END DO
221 :
222 : ! set up matrix S, and needed functions of S
223 882 : CALL ls_scf_init_matrix_s(matrix_s(1)%matrix, ls_scf_env)
224 :
225 : ! get the initial guess for the SCF
226 882 : CALL ls_scf_initial_guess(qs_env, ls_scf_env, nonscf)
227 :
228 882 : IF (ls_scf_env%do_rho_mixing) THEN
229 0 : CALL rho_mixing_ls_init(qs_env, ls_scf_env)
230 : END IF
231 :
232 882 : IF (ls_scf_env%do_pexsi) THEN
233 0 : CALL pexsi_init_scf(ks_env, ls_scf_env%pexsi, matrix_s(1)%matrix)
234 : END IF
235 :
236 882 : IF (qs_env%do_transport) THEN
237 0 : CALL transport_initialize(ks_env, qs_env%transport_env, matrix_s(1)%matrix)
238 : END IF
239 :
240 882 : CALL timestop(handle)
241 :
242 882 : END SUBROUTINE ls_scf_init_scf
243 :
244 : ! **************************************************************************************************
245 : !> \brief deal with the scf initial guess
246 : !> \param qs_env ...
247 : !> \param ls_scf_env ...
248 : !> \param nonscf ...
249 : !> \par History
250 : !> 2012.11 created [Joost VandeVondele]
251 : !> \author Joost VandeVondele
252 : ! **************************************************************************************************
253 1252 : SUBROUTINE ls_scf_initial_guess(qs_env, ls_scf_env, nonscf)
254 : TYPE(qs_environment_type), POINTER :: qs_env
255 : TYPE(ls_scf_env_type) :: ls_scf_env
256 : LOGICAL, INTENT(IN) :: nonscf
257 :
258 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_initial_guess'
259 : INTEGER, PARAMETER :: aspc_guess = 2, atomic_guess = 1, &
260 : restart_guess = 3
261 :
262 : CHARACTER(LEN=default_path_length) :: file_name, project_name
263 : INTEGER :: handle, iaspc, initial_guess_type, &
264 : ispin, istore, naspc, unit_nr
265 : REAL(KIND=dp) :: alpha, cs_pos
266 : TYPE(cp_logger_type), POINTER :: logger
267 : TYPE(dbcsr_distribution_type) :: dist
268 : TYPE(dbcsr_type) :: matrix_tmp1
269 :
270 512 : IF (ls_scf_env%do_pao) RETURN ! pao has its own initial guess
271 :
272 370 : CALL timeset(routineN, handle)
273 :
274 : ! get a useful output_unit
275 370 : logger => cp_get_default_logger()
276 370 : IF (logger%para_env%is_source()) THEN
277 185 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
278 : ELSE
279 : unit_nr = -1
280 : END IF
281 :
282 185 : IF (unit_nr > 0) WRITE (unit_nr, '()')
283 : ! if there is no history go for the atomic guess, otherwise extrapolate the dm history
284 370 : IF (ls_scf_env%scf_history%istore == 0) THEN
285 244 : IF (ls_scf_env%restart_read) THEN
286 : initial_guess_type = restart_guess
287 : ELSE
288 : initial_guess_type = atomic_guess
289 : END IF
290 : ELSE
291 : initial_guess_type = aspc_guess
292 : END IF
293 :
294 : ! how to get the initial guess
295 : SELECT CASE (initial_guess_type)
296 : CASE (atomic_guess)
297 240 : CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env, ls_scf_env%energy_init, nonscf)
298 240 : IF (unit_nr > 0) WRITE (unit_nr, '()')
299 : CASE (restart_guess)
300 4 : project_name = logger%iter_info%project_name
301 8 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
302 4 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_", ispin, "_RESTART.dm"
303 4 : CALL dbcsr_get_info(ls_scf_env%matrix_p(1), distribution=dist)
304 4 : CALL dbcsr_binary_read(file_name, distribution=dist, matrix_new=ls_scf_env%matrix_p(ispin))
305 4 : cs_pos = dbcsr_checksum(ls_scf_env%matrix_p(ispin), pos=.TRUE.)
306 12 : IF (unit_nr > 0) THEN
307 2 : WRITE (unit_nr, '(T2,A,E20.8)') "Read restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
308 : END IF
309 : END DO
310 :
311 : ! directly go to computing the corresponding energy and ks matrix
312 4 : IF (nonscf) THEN
313 0 : CALL ls_nonscf_ks(qs_env, ls_scf_env, ls_scf_env%energy_init)
314 : ELSE
315 4 : CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0)
316 : END IF
317 : CASE (aspc_guess)
318 126 : CALL cite_reference(Kolafa2004)
319 126 : CALL cite_reference(Kuhne2007)
320 126 : naspc = MIN(ls_scf_env%scf_history%istore, ls_scf_env%scf_history%nstore)
321 258 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
322 : ! actual extrapolation
323 132 : CALL dbcsr_set(ls_scf_env%matrix_p(ispin), 0.0_dp)
324 624 : DO iaspc = 1, naspc
325 : alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
326 366 : binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
327 366 : istore = MOD(ls_scf_env%scf_history%istore - iaspc, ls_scf_env%scf_history%nstore) + 1
328 498 : CALL dbcsr_add(ls_scf_env%matrix_p(ispin), ls_scf_env%scf_history%matrix(ispin, istore), 1.0_dp, alpha)
329 : END DO
330 : END DO
331 : END SELECT
332 :
333 : ! which cases need getting purified and non-orthogonal ?
334 126 : SELECT CASE (initial_guess_type)
335 : CASE (atomic_guess, restart_guess)
336 : ! do nothing
337 : CASE (aspc_guess)
338 : ! purification can't be done on the pexsi matrix, which is not necessarily idempotent,
339 : ! and not stored in an ortho basis form
340 126 : IF (.NOT. (ls_scf_env%do_pexsi)) THEN
341 258 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
342 : ! linear combination of P's is not idempotent. A bit of McWeeny is needed to ensure it is again
343 132 : IF (SIZE(ls_scf_env%matrix_p) == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 0.5_dp)
344 : ! to ensure that noisy blocks do not build up during MD (in particular with curvy) filter that guess a bit more
345 132 : CALL dbcsr_filter(ls_scf_env%matrix_p(ispin), ls_scf_env%eps_filter**(2.0_dp/3.0_dp))
346 132 : CALL purify_mcweeny(ls_scf_env%matrix_p(ispin:ispin), ls_scf_env%eps_filter, 3)
347 132 : IF (SIZE(ls_scf_env%matrix_p) == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
348 :
349 258 : IF (ls_scf_env%use_s_sqrt) THEN
350 : ! need to get P in the non-orthogonal basis if it was stored differently
351 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
352 132 : matrix_type=dbcsr_type_no_symmetry)
353 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_p(ispin), &
354 132 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
355 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
356 : 0.0_dp, ls_scf_env%matrix_p(ispin), &
357 132 : filter_eps=ls_scf_env%eps_filter)
358 132 : CALL dbcsr_release(matrix_tmp1)
359 :
360 132 : IF (ls_scf_env%has_s_preconditioner) THEN
361 : CALL apply_matrix_preconditioner(ls_scf_env%matrix_p(ispin), "forward", &
362 126 : ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
363 : END IF
364 : END IF
365 : END DO
366 : END IF
367 :
368 : ! compute corresponding energy and ks matrix
369 496 : IF (nonscf) THEN
370 60 : CALL ls_nonscf_ks(qs_env, ls_scf_env, ls_scf_env%energy_init)
371 : ELSE
372 66 : CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, ls_scf_env%energy_init, iscf=0)
373 : END IF
374 : END SELECT
375 :
376 370 : IF (unit_nr > 0) THEN
377 185 : WRITE (unit_nr, '(T2,A,F20.9)') "Energy with the initial guess:", ls_scf_env%energy_init
378 185 : WRITE (unit_nr, '()')
379 : END IF
380 :
381 370 : CALL timestop(handle)
382 :
383 882 : END SUBROUTINE ls_scf_initial_guess
384 :
385 : ! **************************************************************************************************
386 : !> \brief store a history of matrices for later use in ls_scf_initial_guess
387 : !> \param ls_scf_env ...
388 : !> \par History
389 : !> 2012.11 created [Joost VandeVondele]
390 : !> \author Joost VandeVondele
391 : ! **************************************************************************************************
392 370 : SUBROUTINE ls_scf_store_result(ls_scf_env)
393 : TYPE(ls_scf_env_type) :: ls_scf_env
394 :
395 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_store_result'
396 :
397 : CHARACTER(LEN=default_path_length) :: file_name, project_name
398 : INTEGER :: handle, ispin, istore, unit_nr
399 : REAL(KIND=dp) :: cs_pos
400 : TYPE(cp_logger_type), POINTER :: logger
401 : TYPE(dbcsr_type) :: matrix_tmp1
402 :
403 370 : CALL timeset(routineN, handle)
404 :
405 : ! get a useful output_unit
406 370 : logger => cp_get_default_logger()
407 370 : IF (logger%para_env%is_source()) THEN
408 185 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
409 : ELSE
410 : unit_nr = -1
411 : END IF
412 :
413 370 : IF (ls_scf_env%restart_write) THEN
414 12 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
415 6 : project_name = logger%iter_info%project_name
416 6 : WRITE (file_name, '(A,I0,A)') TRIM(project_name)//"_LS_DM_SPIN_", ispin, "_RESTART.dm"
417 6 : cs_pos = dbcsr_checksum(ls_scf_env%matrix_p(ispin), pos=.TRUE.)
418 6 : IF (unit_nr > 0) THEN
419 3 : WRITE (unit_nr, '(T2,A,E20.8)') "Writing restart DM "//TRIM(file_name)//" with checksum: ", cs_pos
420 : END IF
421 6 : IF (ls_scf_env%do_transport .OR. ls_scf_env%do_pexsi) THEN
422 0 : IF (unit_nr > 0) THEN
423 0 : WRITE (unit_nr, '(T6,A)') "The restart DM "//TRIM(file_name)//" has the sparsity of S, therefore,"
424 0 : WRITE (unit_nr, '(T6,A)') "not compatible with methods that require a full DM! "
425 : END IF
426 : END IF
427 12 : CALL dbcsr_binary_write(ls_scf_env%matrix_p(ispin), file_name)
428 : END DO
429 : END IF
430 :
431 370 : IF (ls_scf_env%scf_history%nstore > 0) THEN
432 362 : ls_scf_env%scf_history%istore = ls_scf_env%scf_history%istore + 1
433 744 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
434 382 : istore = MOD(ls_scf_env%scf_history%istore - 1, ls_scf_env%scf_history%nstore) + 1
435 382 : CALL dbcsr_copy(ls_scf_env%scf_history%matrix(ispin, istore), ls_scf_env%matrix_p(ispin))
436 :
437 : ! if we have the sqrt around, we use it to go to the orthogonal basis
438 744 : IF (ls_scf_env%use_s_sqrt) THEN
439 : ! usually sqrt(S) * P * sqrt(S) should be available, or could be stored at least,
440 : ! so that the next multiplications could be saved.
441 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
442 378 : matrix_type=dbcsr_type_no_symmetry)
443 :
444 378 : IF (ls_scf_env%has_s_preconditioner) THEN
445 : CALL apply_matrix_preconditioner(ls_scf_env%scf_history%matrix(ispin, istore), "backward", &
446 332 : ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
447 : END IF
448 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt, ls_scf_env%scf_history%matrix(ispin, istore), &
449 378 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
450 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt, &
451 : 0.0_dp, ls_scf_env%scf_history%matrix(ispin, istore), &
452 378 : filter_eps=ls_scf_env%eps_filter)
453 378 : CALL dbcsr_release(matrix_tmp1)
454 : END IF
455 :
456 : END DO
457 : END IF
458 :
459 370 : CALL timestop(handle)
460 :
461 370 : END SUBROUTINE ls_scf_store_result
462 :
463 : ! **************************************************************************************************
464 : !> \brief Main SCF routine. Can we keep it clean ?
465 : !> \param qs_env ...
466 : !> \param ls_scf_env ...
467 : !> \param nonscf ...
468 : !> \par History
469 : !> 2010.10 created [Joost VandeVondele]
470 : !> \author Joost VandeVondele
471 : ! **************************************************************************************************
472 882 : SUBROUTINE ls_scf_main(qs_env, ls_scf_env, nonscf)
473 : TYPE(qs_environment_type), POINTER :: qs_env
474 : TYPE(ls_scf_env_type) :: ls_scf_env
475 : LOGICAL, INTENT(IN), OPTIONAL :: nonscf
476 :
477 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_main'
478 :
479 : INTEGER :: handle, iscf, ispin, &
480 : nelectron_spin_real, nmixing, nspin, &
481 : unit_nr
482 : LOGICAL :: check_convergence, diis_step, do_transport, extra_scf, maxscf_reached, &
483 : scf_converged, should_stop, transm_maxscf_reached, transm_scf_converged
484 : REAL(KIND=dp) :: energy_diff, energy_new, energy_old, &
485 : eps_diis, t1, t2, tdiag
486 : TYPE(cp_logger_type), POINTER :: logger
487 882 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
488 882 : TYPE(dbcsr_type), ALLOCATABLE, DIMENSION(:) :: matrix_ks_deviation, matrix_mixing_old
489 : TYPE(energy_correction_type), POINTER :: ec_env
490 : TYPE(qs_diis_buffer_type_sparse), POINTER :: diis_buffer
491 : TYPE(transport_env_type), POINTER :: transport_env
492 :
493 882 : CALL timeset(routineN, handle)
494 :
495 : ! get a useful output_unit
496 882 : logger => cp_get_default_logger()
497 882 : IF (logger%para_env%is_source()) THEN
498 441 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
499 : ELSE
500 441 : unit_nr = -1
501 : END IF
502 :
503 882 : nspin = ls_scf_env%nspins
504 :
505 : ! old quantities, useful for mixing
506 5332 : ALLOCATE (matrix_mixing_old(nspin), matrix_ks_deviation(nspin))
507 1784 : DO ispin = 1, nspin
508 902 : CALL dbcsr_create(matrix_mixing_old(ispin), template=ls_scf_env%matrix_ks(ispin))
509 :
510 902 : CALL dbcsr_create(matrix_ks_deviation(ispin), template=ls_scf_env%matrix_ks(ispin))
511 1784 : CALL dbcsr_set(matrix_ks_deviation(ispin), 0.0_dp)
512 : END DO
513 2646 : ls_scf_env%homo_spin(:) = 0.0_dp
514 2646 : ls_scf_env%lumo_spin(:) = 0.0_dp
515 :
516 882 : transm_scf_converged = .FALSE.
517 882 : transm_maxscf_reached = .FALSE.
518 :
519 882 : energy_old = 0.0_dp
520 882 : IF (ls_scf_env%scf_history%istore > 0) energy_old = ls_scf_env%energy_init
521 882 : check_convergence = .TRUE.
522 882 : iscf = 0
523 882 : IF (ls_scf_env%ls_diis) THEN
524 4 : diis_step = .FALSE.
525 4 : eps_diis = ls_scf_env%eps_diis
526 4 : nmixing = ls_scf_env%nmixing
527 : NULLIFY (diis_buffer)
528 4 : ALLOCATE (diis_buffer)
529 : CALL qs_diis_b_create_sparse(diis_buffer, &
530 4 : nbuffer=ls_scf_env%max_diis)
531 4 : CALL qs_diis_b_clear_sparse(diis_buffer)
532 4 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
533 : END IF
534 :
535 882 : CALL get_qs_env(qs_env, transport_env=transport_env, do_transport=do_transport)
536 :
537 : ! the real SCF loop
538 2878 : DO
539 :
540 : ! check on max SCF or timing/exit
541 2878 : CALL external_control(should_stop, "SCF", start_time=qs_env%start_time, target_time=qs_env%target_time)
542 2878 : IF (do_transport) THEN
543 0 : maxscf_reached = should_stop .OR. iscf >= ls_scf_env%max_scf
544 : ! one extra scf step for post-processing in transmission calculations
545 0 : IF (transport_env%params%method .EQ. transport_transmission) THEN
546 0 : IF (transm_maxscf_reached) THEN
547 0 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') "SCF not converged! "
548 : EXIT
549 : END IF
550 : transm_maxscf_reached = maxscf_reached
551 : ELSE
552 0 : IF (maxscf_reached) THEN
553 0 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') "SCF not converged! "
554 : EXIT
555 : END IF
556 : END IF
557 : ELSE
558 2878 : IF (should_stop .OR. iscf >= ls_scf_env%max_scf) THEN
559 46 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') "SCF not converged! "
560 : ! Skip Harris functional calculation if ground-state is NOT converged
561 46 : IF (qs_env%energy_correction) THEN
562 0 : CALL get_qs_env(qs_env, ec_env=ec_env)
563 0 : IF (ec_env%skip_ec) ec_env%do_skip = .TRUE.
564 : END IF
565 : EXIT
566 : END IF
567 : END IF
568 :
569 2832 : t1 = m_walltime()
570 2832 : iscf = iscf + 1
571 :
572 : ! first get a copy of the current KS matrix
573 2832 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
574 5884 : DO ispin = 1, nspin
575 : CALL matrix_qs_to_ls(ls_scf_env%matrix_ks(ispin), matrix_ks(ispin)%matrix, &
576 3052 : ls_scf_env%ls_mstruct, covariant=.TRUE.)
577 3052 : IF (ls_scf_env%has_s_preconditioner) THEN
578 : CALL apply_matrix_preconditioner(ls_scf_env%matrix_ks(ispin), "forward", &
579 1302 : ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
580 : END IF
581 5884 : CALL dbcsr_filter(ls_scf_env%matrix_ks(ispin), ls_scf_env%eps_filter)
582 : END DO
583 : ! run curvy steps if required. Needs an idempotent DM (either perification or restart)
584 2832 : IF ((iscf > 1 .OR. ls_scf_env%scf_history%istore > 0) .AND. ls_scf_env%curvy_steps) THEN
585 90 : CALL dm_ls_curvy_optimization(ls_scf_env, energy_old, check_convergence)
586 : ELSE
587 : ! turn the KS matrix in a density matrix
588 5692 : DO ispin = 1, nspin
589 2950 : IF (nonscf) THEN
590 66 : CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin))
591 2884 : ELSEIF (ls_scf_env%do_rho_mixing) THEN
592 0 : CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin))
593 : ELSE
594 2884 : IF (iscf == 1) THEN
595 : ! initialize the mixing matrix with the current state if needed
596 826 : CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin))
597 : ELSE
598 2058 : IF (ls_scf_env%ls_diis) THEN ! ------- IF-DIIS+MIX--- START
599 28 : IF (diis_step .AND. (iscf - 1) .GE. ls_scf_env%iter_ini_diis) THEN
600 22 : IF (unit_nr > 0) THEN
601 : WRITE (unit_nr, '(A61)') &
602 11 : '*************************************************************'
603 : WRITE (unit_nr, '(A50,2(I3,A1),L1,A1)') &
604 11 : " Using DIIS mixed KS: (iscf,INI_DIIS,DIIS_STEP)=(", &
605 22 : iscf, ",", ls_scf_env%iter_ini_diis, ",", diis_step, ")"
606 : WRITE (unit_nr, '(A52)') &
607 11 : " KS_nw= DIIS-Linear-Combination-Previous KS matrices"
608 : WRITE (unit_nr, '(61A)') &
609 11 : "*************************************************************"
610 : END IF
611 : CALL dbcsr_copy(matrix_mixing_old(ispin), & ! out
612 22 : ls_scf_env%matrix_ks(ispin)) ! in
613 : ELSE
614 6 : IF (unit_nr > 0) THEN
615 : WRITE (unit_nr, '(A57)') &
616 3 : "*********************************************************"
617 : WRITE (unit_nr, '(A23,F5.3,A25,I3)') &
618 3 : " Using MIXING_FRACTION=", ls_scf_env%mixing_fraction, &
619 6 : " to mix KS matrix: iscf=", iscf
620 : WRITE (unit_nr, '(A7,F5.3,A6,F5.3,A7)') &
621 3 : " KS_nw=", ls_scf_env%mixing_fraction, "*KS + ", &
622 6 : 1.0_dp - ls_scf_env%mixing_fraction, "*KS_old"
623 : WRITE (unit_nr, '(A57)') &
624 3 : "*********************************************************"
625 : END IF
626 : ! perform the mixing of ks matrices
627 : CALL dbcsr_add(matrix_mixing_old(ispin), &
628 : ls_scf_env%matrix_ks(ispin), &
629 : 1.0_dp - ls_scf_env%mixing_fraction, &
630 6 : ls_scf_env%mixing_fraction)
631 : END IF
632 : ELSE ! otherwise
633 2030 : IF (unit_nr > 0) THEN
634 : WRITE (unit_nr, '(A57)') &
635 1015 : "*********************************************************"
636 : WRITE (unit_nr, '(A23,F5.3,A25,I3)') &
637 1015 : " Using MIXING_FRACTION=", ls_scf_env%mixing_fraction, &
638 2030 : " to mix KS matrix: iscf=", iscf
639 : WRITE (unit_nr, '(A7,F5.3,A6,F5.3,A7)') &
640 1015 : " KS_nw=", ls_scf_env%mixing_fraction, "*KS + ", &
641 2030 : 1.0_dp - ls_scf_env%mixing_fraction, "*KS_old"
642 : WRITE (unit_nr, '(A57)') &
643 1015 : "*********************************************************"
644 : END IF
645 : ! perform the mixing of ks matrices
646 : CALL dbcsr_add(matrix_mixing_old(ispin), &
647 : ls_scf_env%matrix_ks(ispin), &
648 : 1.0_dp - ls_scf_env%mixing_fraction, &
649 2030 : ls_scf_env%mixing_fraction)
650 : END IF ! ------- IF-DIIS+MIX--- END
651 : END IF
652 : END IF
653 :
654 : ! compute the density matrix that matches it
655 : ! we need the proper number of states
656 2950 : nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
657 2950 : IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
658 :
659 2950 : IF (do_transport) THEN
660 0 : IF (ls_scf_env%has_s_preconditioner) &
661 0 : CPABORT("NOT YET IMPLEMENTED with S preconditioner. ")
662 0 : IF (ls_scf_env%ls_mstruct%cluster_type .NE. ls_cluster_atomic) &
663 0 : CPABORT("NOT YET IMPLEMENTED with molecular clustering. ")
664 :
665 0 : extra_scf = maxscf_reached .OR. scf_converged
666 : ! get the current Kohn-Sham matrix (ks) and return matrix_p evaluated using an external C routine
667 : CALL external_scf_method(transport_env, ls_scf_env%matrix_s, matrix_mixing_old(ispin), &
668 : ls_scf_env%matrix_p(ispin), nelectron_spin_real, ls_scf_env%natoms, &
669 0 : energy_diff, iscf, extra_scf)
670 :
671 : ELSE
672 3858 : SELECT CASE (ls_scf_env%purification_method)
673 : CASE (ls_scf_sign)
674 : CALL density_matrix_sign(ls_scf_env%matrix_p(ispin), ls_scf_env%mu_spin(ispin), ls_scf_env%fixed_mu, &
675 : ls_scf_env%sign_method, ls_scf_env%sign_order, matrix_mixing_old(ispin), &
676 : ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, nelectron_spin_real, &
677 : ls_scf_env%eps_filter, ls_scf_env%sign_symmetric, ls_scf_env%submatrix_sign_method, &
678 908 : ls_scf_env%matrix_s_sqrt_inv)
679 : CASE (ls_scf_tc2)
680 : CALL density_matrix_tc2(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, &
681 : nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), &
682 : ls_scf_env%lumo_spin(ispin), non_monotonic=ls_scf_env%non_monotonic, &
683 : eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
684 118 : iounit=-1)
685 : CASE (ls_scf_trs4)
686 : CALL density_matrix_trs4(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, &
687 : nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), &
688 : ls_scf_env%lumo_spin(ispin), ls_scf_env%mu_spin(ispin), &
689 : dynamic_threshold=ls_scf_env%dynamic_threshold, &
690 : matrix_ks_deviation=matrix_ks_deviation(ispin), &
691 : eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
692 1924 : iounit=-1)
693 : CASE (ls_scf_pexsi)
694 0 : IF (ls_scf_env%has_s_preconditioner) &
695 0 : CPABORT("S preconditioning not implemented in combination with the PEXSI library. ")
696 0 : IF (ls_scf_env%ls_mstruct%cluster_type .NE. ls_cluster_atomic) &
697 : CALL cp_abort(__LOCATION__, &
698 0 : "Molecular clustering not implemented in combination with the PEXSI library. ")
699 : CALL density_matrix_pexsi(ls_scf_env%pexsi, ls_scf_env%matrix_p(ispin), ls_scf_env%pexsi%matrix_w(ispin), &
700 : ls_scf_env%pexsi%kTS(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s, &
701 2950 : nelectron_spin_real, ls_scf_env%mu_spin(ispin), iscf, ispin)
702 : END SELECT
703 : END IF
704 :
705 2950 : IF (ls_scf_env%has_s_preconditioner) THEN
706 : CALL apply_matrix_preconditioner(ls_scf_env%matrix_p(ispin), "forward", &
707 1302 : ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
708 : END IF
709 2950 : CALL dbcsr_filter(ls_scf_env%matrix_p(ispin), ls_scf_env%eps_filter)
710 :
711 5692 : IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
712 :
713 : END DO
714 : END IF
715 :
716 : ! compute the corresponding new energy KS matrix and new energy
717 2832 : IF (nonscf) THEN
718 66 : CALL ls_nonscf_energy(qs_env, ls_scf_env)
719 : ELSE
720 2766 : CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf)
721 : END IF
722 :
723 2832 : IF (ls_scf_env%do_pexsi) THEN
724 0 : CALL pexsi_to_qs(ls_scf_env, qs_env, kTS=ls_scf_env%pexsi%kTS)
725 : END IF
726 :
727 2832 : t2 = m_walltime()
728 2832 : IF (nonscf) THEN
729 66 : tdiag = t2 - t1
730 66 : CALL qs_nonscf_print_summary(qs_env, tdiag, ls_scf_env%nelectron_total, unit_nr)
731 66 : EXIT
732 : ELSE
733 : ! report current SCF loop
734 2766 : energy_diff = energy_new - energy_old
735 2766 : energy_old = energy_new
736 2766 : IF (unit_nr > 0) THEN
737 1383 : WRITE (unit_nr, *)
738 1383 : WRITE (unit_nr, '(T2,A,I6,F20.9,F20.9,F12.6)') "SCF", iscf, energy_new, energy_diff, t2 - t1
739 1383 : WRITE (unit_nr, *)
740 1383 : CALL m_flush(unit_nr)
741 : END IF
742 : END IF
743 :
744 2766 : IF (do_transport) THEN
745 0 : scf_converged = check_convergence .AND. ABS(energy_diff) < ls_scf_env%eps_scf*ls_scf_env%nelectron_total
746 : ! one extra scf step for post-processing in transmission calculations
747 0 : IF (transport_env%params%method .EQ. transport_transmission) THEN
748 0 : IF (transm_scf_converged) EXIT
749 : transm_scf_converged = scf_converged
750 : ELSE
751 0 : IF (scf_converged) THEN
752 0 : IF (unit_nr > 0) WRITE (unit_nr, '(/,T2,A,I5,A/)') "SCF run converged in ", iscf, " steps."
753 : EXIT
754 : END IF
755 : END IF
756 : ELSE
757 : ! exit criterion on the energy only for the time being
758 2766 : IF (check_convergence .AND. ABS(energy_diff) < ls_scf_env%eps_scf*ls_scf_env%nelectron_total) THEN
759 770 : IF (unit_nr > 0) WRITE (unit_nr, '(/,T2,A,I5,A/)') "SCF run converged in ", iscf, " steps."
760 : ! Skip Harris functional calculation if ground-state is NOT converged
761 770 : IF (qs_env%energy_correction) THEN
762 20 : CALL get_qs_env(qs_env, ec_env=ec_env)
763 20 : IF (ec_env%skip_ec) ec_env%do_skip = .FALSE.
764 : END IF
765 : EXIT
766 : END IF
767 : END IF
768 :
769 1996 : IF (ls_scf_env%ls_diis) THEN
770 : ! diis_buffer, buffer with 1) Kohn-Sham history matrix,
771 : ! 2) KS error history matrix (f=KPS-SPK),
772 : ! 3) B matrix (for finding DIIS weighting coefficients)
773 : CALL qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, &
774 : iscf, diis_step, eps_diis, nmixing, matrix_s(1)%matrix, &
775 18 : ls_scf_env%eps_filter)
776 : END IF
777 :
778 1996 : IF (ls_scf_env%do_pexsi) THEN
779 : CALL pexsi_set_convergence_tolerance(ls_scf_env%pexsi, energy_diff, &
780 : ls_scf_env%eps_scf*ls_scf_env%nelectron_total, &
781 : ! initialize in second scf step of first SCF cycle:
782 : (iscf .EQ. 2) .AND. (ls_scf_env%scf_history%istore .EQ. 0), &
783 0 : check_convergence)
784 : END IF
785 :
786 : END DO
787 :
788 : ! free storage
789 882 : IF (ls_scf_env%ls_diis) THEN
790 4 : CALL qs_diis_b_release_sparse(diis_buffer)
791 4 : DEALLOCATE (diis_buffer)
792 : END IF
793 1784 : DO ispin = 1, nspin
794 902 : CALL dbcsr_release(matrix_mixing_old(ispin))
795 1784 : CALL dbcsr_release(matrix_ks_deviation(ispin))
796 : END DO
797 882 : DEALLOCATE (matrix_mixing_old, matrix_ks_deviation)
798 :
799 882 : CALL timestop(handle)
800 :
801 882 : END SUBROUTINE ls_scf_main
802 :
803 : ! **************************************************************************************************
804 : !> \brief after SCF we have a density matrix, and the self consistent KS matrix
805 : !> analyze its properties.
806 : !> \param qs_env ...
807 : !> \param ls_scf_env ...
808 : !> \par History
809 : !> 2010.10 created [Joost VandeVondele]
810 : !> \author Joost VandeVondele
811 : ! **************************************************************************************************
812 882 : SUBROUTINE ls_scf_post(qs_env, ls_scf_env)
813 : TYPE(qs_environment_type), POINTER :: qs_env
814 : TYPE(ls_scf_env_type) :: ls_scf_env
815 :
816 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_post'
817 :
818 : INTEGER :: handle, ispin, unit_nr
819 : REAL(KIND=dp) :: occ
820 : TYPE(cp_logger_type), POINTER :: logger
821 882 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
822 : TYPE(dft_control_type), POINTER :: dft_control
823 :
824 882 : CALL timeset(routineN, handle)
825 :
826 882 : CALL get_qs_env(qs_env, dft_control=dft_control)
827 :
828 : ! get a useful output_unit
829 882 : logger => cp_get_default_logger()
830 882 : IF (logger%para_env%is_source()) THEN
831 441 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
832 : ELSE
833 441 : unit_nr = -1
834 : END IF
835 :
836 : ! store the matrix for a next scf run
837 882 : IF (.NOT. ls_scf_env%do_pao) &
838 370 : CALL ls_scf_store_result(ls_scf_env)
839 :
840 : ! write homo and lumo energy and occupation (if not already part of the output)
841 882 : IF (ls_scf_env%curvy_steps) THEN
842 18 : CALL post_scf_homo_lumo(ls_scf_env)
843 :
844 : ! always report P occ
845 18 : IF (unit_nr > 0) WRITE (unit_nr, *) ""
846 38 : DO ispin = 1, ls_scf_env%nspins
847 20 : occ = dbcsr_get_occupation(ls_scf_env%matrix_p(ispin))
848 38 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Density matrix (P) occupation ", occ
849 : END DO
850 : END IF
851 :
852 : ! compute the matrix_w if associated
853 882 : IF (ls_scf_env%calculate_forces) THEN
854 170 : CALL get_qs_env(qs_env, matrix_w=matrix_w)
855 170 : CPASSERT(ASSOCIATED(matrix_w))
856 170 : IF (ls_scf_env%do_pexsi) THEN
857 0 : CALL pexsi_to_qs(ls_scf_env, qs_env, matrix_w=ls_scf_env%pexsi%matrix_w)
858 : ELSE
859 170 : CALL calculate_w_matrix_ls(matrix_w, ls_scf_env)
860 : END IF
861 : END IF
862 :
863 : ! compute properties
864 :
865 882 : IF (ls_scf_env%perform_mu_scan) CALL post_scf_mu_scan(ls_scf_env)
866 :
867 882 : IF (ls_scf_env%report_all_sparsities) CALL post_scf_sparsities(ls_scf_env)
868 :
869 882 : IF (dft_control%qs_control%dftb) THEN
870 44 : CALL scf_post_calculation_tb(qs_env, "DFTB", .TRUE.)
871 838 : ELSE IF (dft_control%qs_control%xtb) THEN
872 86 : CALL scf_post_calculation_tb(qs_env, "xTB", .TRUE.)
873 : ELSE
874 752 : CALL write_mo_free_results(qs_env)
875 : END IF
876 :
877 882 : IF (ls_scf_env%chebyshev%compute_chebyshev) CALL compute_chebyshev(qs_env, ls_scf_env)
878 :
879 882 : IF (.TRUE.) CALL post_scf_experiment()
880 :
881 882 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
882 : !
883 : ELSE
884 752 : CALL qs_scf_post_moments(qs_env%input, logger, qs_env, unit_nr)
885 : END IF
886 :
887 : ! clean up used data
888 :
889 882 : CALL dbcsr_release(ls_scf_env%matrix_s)
890 882 : CALL deallocate_curvy_data(ls_scf_env%curvy_data)
891 :
892 882 : IF (ls_scf_env%has_s_preconditioner) THEN
893 326 : CALL dbcsr_release(ls_scf_env%matrix_bs_sqrt)
894 326 : CALL dbcsr_release(ls_scf_env%matrix_bs_sqrt_inv)
895 : END IF
896 :
897 882 : IF (ls_scf_env%needs_s_inv) THEN
898 880 : CALL dbcsr_release(ls_scf_env%matrix_s_inv)
899 : END IF
900 :
901 882 : IF (ls_scf_env%use_s_sqrt) THEN
902 878 : CALL dbcsr_release(ls_scf_env%matrix_s_sqrt)
903 878 : CALL dbcsr_release(ls_scf_env%matrix_s_sqrt_inv)
904 : END IF
905 :
906 1784 : DO ispin = 1, SIZE(ls_scf_env%matrix_ks)
907 1784 : CALL dbcsr_release(ls_scf_env%matrix_ks(ispin))
908 : END DO
909 882 : DEALLOCATE (ls_scf_env%matrix_ks)
910 :
911 882 : IF (ls_scf_env%do_pexsi) &
912 0 : CALL pexsi_finalize_scf(ls_scf_env%pexsi, ls_scf_env%mu_spin)
913 :
914 882 : CALL timestop(handle)
915 :
916 882 : END SUBROUTINE ls_scf_post
917 :
918 : ! **************************************************************************************************
919 : !> \brief Compute the HOMO LUMO energies post SCF
920 : !> \param ls_scf_env ...
921 : !> \par History
922 : !> 2013.06 created [Joost VandeVondele]
923 : !> \author Joost VandeVondele
924 : ! **************************************************************************************************
925 18 : SUBROUTINE post_scf_homo_lumo(ls_scf_env)
926 : TYPE(ls_scf_env_type) :: ls_scf_env
927 :
928 : CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_homo_lumo'
929 :
930 : INTEGER :: handle, ispin, nspin, unit_nr
931 : LOGICAL :: converged
932 : REAL(KIND=dp) :: eps_max, eps_min, homo, lumo
933 : TYPE(cp_logger_type), POINTER :: logger
934 : TYPE(dbcsr_type) :: matrix_k, matrix_p, matrix_tmp
935 :
936 18 : CALL timeset(routineN, handle)
937 :
938 : ! get a useful output_unit
939 18 : logger => cp_get_default_logger()
940 18 : IF (logger%para_env%is_source()) THEN
941 9 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
942 : ELSE
943 9 : unit_nr = -1
944 : END IF
945 :
946 18 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') ""
947 :
948 : ! TODO: remove these limitations
949 18 : CPASSERT(.NOT. ls_scf_env%has_s_preconditioner)
950 18 : CPASSERT(ls_scf_env%use_s_sqrt)
951 :
952 18 : nspin = ls_scf_env%nspins
953 :
954 18 : CALL dbcsr_create(matrix_p, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
955 :
956 18 : CALL dbcsr_create(matrix_k, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
957 :
958 18 : CALL dbcsr_create(matrix_tmp, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
959 :
960 38 : DO ispin = 1, nspin
961 : ! ortho basis ks
962 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_ks(ispin), &
963 20 : 0.0_dp, matrix_tmp, filter_eps=ls_scf_env%eps_filter)
964 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, ls_scf_env%matrix_s_sqrt_inv, &
965 20 : 0.0_dp, matrix_k, filter_eps=ls_scf_env%eps_filter)
966 :
967 : ! extremal eigenvalues ks
968 : CALL arnoldi_extremal(matrix_k, eps_max, eps_min, max_iter=ls_scf_env%max_iter_lanczos, &
969 20 : threshold=ls_scf_env%eps_lanczos, converged=converged)
970 :
971 : ! ortho basis p
972 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_p(ispin), &
973 20 : 0.0_dp, matrix_tmp, filter_eps=ls_scf_env%eps_filter)
974 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, ls_scf_env%matrix_s_sqrt, &
975 20 : 0.0_dp, matrix_p, filter_eps=ls_scf_env%eps_filter)
976 20 : IF (nspin == 1) CALL dbcsr_scale(matrix_p, 0.5_dp)
977 :
978 : ! go compute homo lumo
979 : CALL compute_homo_lumo(matrix_k, matrix_p, eps_min, eps_max, ls_scf_env%eps_filter, &
980 58 : ls_scf_env%max_iter_lanczos, ls_scf_env%eps_lanczos, homo, lumo, unit_nr)
981 :
982 : END DO
983 :
984 18 : CALL dbcsr_release(matrix_p)
985 18 : CALL dbcsr_release(matrix_k)
986 18 : CALL dbcsr_release(matrix_tmp)
987 :
988 18 : CALL timestop(handle)
989 :
990 18 : END SUBROUTINE post_scf_homo_lumo
991 :
992 : ! **************************************************************************************************
993 : !> \brief Compute the density matrix for various values of the chemical potential
994 : !> \param ls_scf_env ...
995 : !> \par History
996 : !> 2010.10 created [Joost VandeVondele]
997 : !> \author Joost VandeVondele
998 : ! **************************************************************************************************
999 2 : SUBROUTINE post_scf_mu_scan(ls_scf_env)
1000 : TYPE(ls_scf_env_type) :: ls_scf_env
1001 :
1002 : CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_mu_scan'
1003 :
1004 : INTEGER :: handle, imu, ispin, nelectron_spin_real, &
1005 : nmu, nspin, unit_nr
1006 : REAL(KIND=dp) :: mu, t1, t2, trace
1007 : TYPE(cp_logger_type), POINTER :: logger
1008 : TYPE(dbcsr_type) :: matrix_p
1009 :
1010 2 : CALL timeset(routineN, handle)
1011 :
1012 : ! get a useful output_unit
1013 2 : logger => cp_get_default_logger()
1014 2 : IF (logger%para_env%is_source()) THEN
1015 1 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1016 : ELSE
1017 : unit_nr = -1
1018 : END IF
1019 :
1020 2 : nspin = ls_scf_env%nspins
1021 :
1022 2 : CALL dbcsr_create(matrix_p, template=ls_scf_env%matrix_p(1))
1023 :
1024 2 : nmu = 10
1025 24 : DO imu = 0, nmu
1026 :
1027 22 : t1 = m_walltime()
1028 :
1029 22 : mu = -0.4_dp + imu*(0.1_dp + 0.4_dp)/nmu
1030 :
1031 22 : IF (unit_nr > 0) WRITE (unit_nr, *) "------- starting with mu ", mu
1032 :
1033 44 : DO ispin = 1, nspin
1034 : ! we need the proper number of states
1035 22 : nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
1036 22 : IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
1037 :
1038 : CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, ls_scf_env%sign_method, &
1039 : ls_scf_env%sign_order, ls_scf_env%matrix_ks(ispin), &
1040 : ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, &
1041 : ls_scf_env%eps_filter, ls_scf_env%sign_symmetric, &
1042 44 : ls_scf_env%submatrix_sign_method, ls_scf_env%matrix_s_sqrt_inv)
1043 : END DO
1044 :
1045 22 : t2 = m_walltime()
1046 :
1047 24 : IF (unit_nr > 0) WRITE (unit_nr, *) " obtained ", mu, trace, t2 - t1
1048 :
1049 : END DO
1050 :
1051 2 : CALL dbcsr_release(matrix_p)
1052 :
1053 2 : CALL timestop(handle)
1054 :
1055 2 : END SUBROUTINE post_scf_mu_scan
1056 :
1057 : ! **************************************************************************************************
1058 : !> \brief Report on the sparsities of various interesting matrices.
1059 : !>
1060 : !> \param ls_scf_env ...
1061 : !> \par History
1062 : !> 2010.10 created [Joost VandeVondele]
1063 : !> \author Joost VandeVondele
1064 : ! **************************************************************************************************
1065 180 : SUBROUTINE post_scf_sparsities(ls_scf_env)
1066 : TYPE(ls_scf_env_type) :: ls_scf_env
1067 :
1068 : CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_sparsities'
1069 :
1070 : CHARACTER(LEN=default_string_length) :: title
1071 : INTEGER :: handle, ispin, nspin, unit_nr
1072 : TYPE(cp_logger_type), POINTER :: logger
1073 : TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2
1074 :
1075 180 : CALL timeset(routineN, handle)
1076 :
1077 : ! get a useful output_unit
1078 180 : logger => cp_get_default_logger()
1079 180 : IF (logger%para_env%is_source()) THEN
1080 90 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1081 : ELSE
1082 90 : unit_nr = -1
1083 : END IF
1084 :
1085 180 : nspin = ls_scf_env%nspins
1086 :
1087 180 : IF (unit_nr > 0) THEN
1088 90 : WRITE (unit_nr, '()')
1089 90 : WRITE (unit_nr, '(T2,A,E17.3)') "Sparsity reports for eps_filter: ", ls_scf_env%eps_filter
1090 90 : WRITE (unit_nr, '()')
1091 : END IF
1092 :
1093 : CALL report_matrix_sparsity(ls_scf_env%matrix_s, unit_nr, "overlap matrix (S)", &
1094 180 : ls_scf_env%eps_filter)
1095 :
1096 368 : DO ispin = 1, nspin
1097 188 : WRITE (title, '(A,I3)') "Kohn-Sham matrix (H) for spin ", ispin
1098 : CALL report_matrix_sparsity(ls_scf_env%matrix_ks(ispin), unit_nr, title, &
1099 368 : ls_scf_env%eps_filter)
1100 : END DO
1101 :
1102 180 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1103 180 : CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1104 :
1105 368 : DO ispin = 1, nspin
1106 188 : WRITE (title, '(A,I3)') "Density matrix (P) for spin ", ispin
1107 : CALL report_matrix_sparsity(ls_scf_env%matrix_p(ispin), unit_nr, title, &
1108 188 : ls_scf_env%eps_filter)
1109 :
1110 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s, ls_scf_env%matrix_p(ispin), &
1111 188 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1112 :
1113 188 : WRITE (title, '(A,I3,A)') "S * P(", ispin, ")"
1114 188 : CALL report_matrix_sparsity(matrix_tmp1, unit_nr, title, ls_scf_env%eps_filter)
1115 :
1116 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s, &
1117 188 : 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
1118 188 : WRITE (title, '(A,I3,A)') "S * P(", ispin, ") * S"
1119 368 : CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
1120 : END DO
1121 :
1122 180 : IF (ls_scf_env%needs_s_inv) THEN
1123 : CALL report_matrix_sparsity(ls_scf_env%matrix_s_inv, unit_nr, "inv(S)", &
1124 180 : ls_scf_env%eps_filter)
1125 368 : DO ispin = 1, nspin
1126 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_inv, ls_scf_env%matrix_ks(ispin), &
1127 188 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1128 :
1129 188 : WRITE (title, '(A,I3,A)') "inv(S) * H(", ispin, ")"
1130 368 : CALL report_matrix_sparsity(matrix_tmp1, unit_nr, title, ls_scf_env%eps_filter)
1131 : END DO
1132 : END IF
1133 :
1134 180 : IF (ls_scf_env%use_s_sqrt) THEN
1135 :
1136 : CALL report_matrix_sparsity(ls_scf_env%matrix_s_sqrt, unit_nr, "sqrt(S)", &
1137 178 : ls_scf_env%eps_filter)
1138 : CALL report_matrix_sparsity(ls_scf_env%matrix_s_sqrt_inv, unit_nr, "inv(sqrt(S))", &
1139 178 : ls_scf_env%eps_filter)
1140 :
1141 364 : DO ispin = 1, nspin
1142 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_ks(ispin), &
1143 186 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1144 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
1145 186 : 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
1146 186 : WRITE (title, '(A,I3,A)') "inv(sqrt(S)) * H(", ispin, ") * inv(sqrt(S))"
1147 364 : CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
1148 : END DO
1149 :
1150 364 : DO ispin = 1, nspin
1151 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_p(ispin), &
1152 186 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1153 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt, &
1154 186 : 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
1155 186 : WRITE (title, '(A,I3,A)') "sqrt(S) * P(", ispin, ") * sqrt(S)"
1156 364 : CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
1157 : END DO
1158 :
1159 : END IF
1160 :
1161 180 : CALL dbcsr_release(matrix_tmp1)
1162 180 : CALL dbcsr_release(matrix_tmp2)
1163 :
1164 180 : CALL timestop(handle)
1165 :
1166 180 : END SUBROUTINE post_scf_sparsities
1167 :
1168 : ! **************************************************************************************************
1169 : !> \brief Helper routine to report on the sparsity of a single matrix,
1170 : !> for several filtering values
1171 : !> \param matrix ...
1172 : !> \param unit_nr ...
1173 : !> \param title ...
1174 : !> \param eps ...
1175 : !> \par History
1176 : !> 2010.10 created [Joost VandeVondele]
1177 : !> \author Joost VandeVondele
1178 : ! **************************************************************************************************
1179 2028 : SUBROUTINE report_matrix_sparsity(matrix, unit_nr, title, eps)
1180 : TYPE(dbcsr_type) :: matrix
1181 : INTEGER :: unit_nr
1182 : CHARACTER(LEN=*) :: title
1183 : REAL(KIND=dp) :: eps
1184 :
1185 : CHARACTER(len=*), PARAMETER :: routineN = 'report_matrix_sparsity'
1186 :
1187 : INTEGER :: handle
1188 : REAL(KIND=dp) :: eps_local, occ
1189 : TYPE(dbcsr_type) :: matrix_tmp
1190 :
1191 2028 : CALL timeset(routineN, handle)
1192 2028 : CALL dbcsr_create(matrix_tmp, template=matrix, name=TRIM(title))
1193 2028 : CALL dbcsr_copy(matrix_tmp, matrix, name=TRIM(title))
1194 :
1195 2028 : IF (unit_nr > 0) THEN
1196 1014 : WRITE (unit_nr, '(T2,A)') "Sparsity for : "//TRIM(title)
1197 : END IF
1198 :
1199 2028 : eps_local = MAX(eps, 10E-14_dp)
1200 14796 : DO
1201 16824 : IF (eps_local > 1.1_dp) EXIT
1202 14796 : CALL dbcsr_filter(matrix_tmp, eps_local)
1203 14796 : occ = dbcsr_get_occupation(matrix_tmp)
1204 14796 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,F16.12,A3,F16.12)') eps_local, " : ", occ
1205 14796 : eps_local = eps_local*10
1206 : END DO
1207 :
1208 2028 : CALL dbcsr_release(matrix_tmp)
1209 :
1210 2028 : CALL timestop(handle)
1211 :
1212 2028 : END SUBROUTINE report_matrix_sparsity
1213 :
1214 : ! **************************************************************************************************
1215 : !> \brief Compute matrix_w as needed for the forces
1216 : !> \param matrix_w ...
1217 : !> \param ls_scf_env ...
1218 : !> \par History
1219 : !> 2010.11 created [Joost VandeVondele]
1220 : !> \author Joost VandeVondele
1221 : ! **************************************************************************************************
1222 200 : SUBROUTINE calculate_w_matrix_ls(matrix_w, ls_scf_env)
1223 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
1224 : TYPE(ls_scf_env_type) :: ls_scf_env
1225 :
1226 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_ls'
1227 :
1228 : INTEGER :: handle, ispin
1229 : REAL(KIND=dp) :: scaling
1230 : TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2, matrix_tmp3
1231 :
1232 200 : CALL timeset(routineN, handle)
1233 :
1234 200 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1235 200 : CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1236 200 : CALL dbcsr_create(matrix_tmp3, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1237 :
1238 200 : IF (ls_scf_env%nspins == 1) THEN
1239 192 : scaling = 0.5_dp
1240 : ELSE
1241 8 : scaling = 1.0_dp
1242 : END IF
1243 :
1244 408 : DO ispin = 1, ls_scf_env%nspins
1245 :
1246 208 : CALL dbcsr_copy(matrix_tmp3, ls_scf_env%matrix_ks(ispin))
1247 208 : IF (ls_scf_env%has_s_preconditioner) THEN
1248 : CALL apply_matrix_preconditioner(matrix_tmp3, "backward", &
1249 136 : ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
1250 : END IF
1251 208 : CALL dbcsr_filter(matrix_tmp3, ls_scf_env%eps_filter)
1252 :
1253 : CALL dbcsr_multiply("N", "N", scaling, ls_scf_env%matrix_p(ispin), matrix_tmp3, &
1254 208 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1255 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_p(ispin), &
1256 208 : 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
1257 408 : CALL matrix_ls_to_qs(matrix_w(ispin)%matrix, matrix_tmp2, ls_scf_env%ls_mstruct, covariant=.FALSE.)
1258 : END DO
1259 :
1260 200 : CALL dbcsr_release(matrix_tmp1)
1261 200 : CALL dbcsr_release(matrix_tmp2)
1262 200 : CALL dbcsr_release(matrix_tmp3)
1263 :
1264 200 : CALL timestop(handle)
1265 :
1266 200 : END SUBROUTINE calculate_w_matrix_ls
1267 :
1268 : ! **************************************************************************************************
1269 : !> \brief a place for quick experiments
1270 : !> \par History
1271 : !> 2010.11 created [Joost VandeVondele]
1272 : !> \author Joost VandeVondele
1273 : ! **************************************************************************************************
1274 882 : SUBROUTINE post_scf_experiment()
1275 :
1276 : CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_experiment'
1277 :
1278 : INTEGER :: handle, unit_nr
1279 : TYPE(cp_logger_type), POINTER :: logger
1280 :
1281 882 : CALL timeset(routineN, handle)
1282 :
1283 : ! get a useful output_unit
1284 882 : logger => cp_get_default_logger()
1285 882 : IF (logger%para_env%is_source()) THEN
1286 441 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1287 : ELSE
1288 : unit_nr = -1
1289 : END IF
1290 :
1291 882 : CALL timestop(handle)
1292 882 : END SUBROUTINE post_scf_experiment
1293 :
1294 : END MODULE dm_ls_scf
|
