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 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 764 : 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 764 : CALL timeset(routineN, handle)
118 764 : do_scf = .TRUE.
119 764 : IF (PRESENT(nonscf)) do_scf = .NOT. nonscf
120 :
121 : ! Moved here from qs_environment to remove dependencies
122 764 : CALL ls_scf_create(qs_env)
123 764 : CALL get_qs_env(qs_env, ls_scf_env=ls_scf_env)
124 :
125 764 : IF (do_scf) THEN
126 698 : CALL pao_optimization_start(qs_env, ls_scf_env)
127 698 : pao_is_done = .FALSE.
128 1614 : DO WHILE (.NOT. pao_is_done)
129 916 : CALL ls_scf_init_scf(qs_env, ls_scf_env, .FALSE.)
130 916 : CALL pao_update(qs_env, ls_scf_env, pao_is_done)
131 916 : CALL ls_scf_main(qs_env, ls_scf_env, .FALSE.)
132 916 : CALL pao_post_scf(qs_env, ls_scf_env, pao_is_done)
133 916 : CALL ls_scf_post(qs_env, ls_scf_env)
134 : END DO
135 698 : 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 764 : CALL timestop(handle)
143 :
144 764 : 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 982 : 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 982 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, matrix_w
165 : TYPE(dft_control_type), POINTER :: dft_control
166 982 : 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 982 : CALL timeset(routineN, handle)
171 :
172 : ! get a useful output_unit
173 982 : logger => cp_get_default_logger()
174 982 : IF (logger%para_env%is_source()) THEN
175 491 : 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 982 : nelectron_spin=ls_scf_env%nelectron_spin)
191 :
192 : ! needs forces ? There might be a better way to flag this
193 982 : ls_scf_env%calculate_forces = ASSOCIATED(matrix_w)
194 :
195 : ! some basic initialization of the QS side of things
196 982 : 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 982 : ls_mstruct=ls_scf_env%ls_mstruct)
201 :
202 982 : nspin = ls_scf_env%nspins
203 982 : IF (ALLOCATED(ls_scf_env%matrix_p)) THEN
204 1198 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
205 1198 : CALL dbcsr_release(ls_scf_env%matrix_p(ispin))
206 : END DO
207 : ELSE
208 1566 : ALLOCATE (ls_scf_env%matrix_p(nspin))
209 : END IF
210 :
211 1992 : DO ispin = 1, nspin
212 : CALL dbcsr_create(ls_scf_env%matrix_p(ispin), template=ls_scf_env%matrix_s, &
213 1992 : matrix_type=dbcsr_type_no_symmetry)
214 : END DO
215 :
216 3956 : ALLOCATE (ls_scf_env%matrix_ks(nspin))
217 1992 : DO ispin = 1, nspin
218 : CALL dbcsr_create(ls_scf_env%matrix_ks(ispin), template=ls_scf_env%matrix_s, &
219 1992 : matrix_type=dbcsr_type_no_symmetry)
220 : END DO
221 :
222 : ! set up matrix S, and needed functions of S
223 982 : CALL ls_scf_init_matrix_s(matrix_s(1)%matrix, ls_scf_env)
224 :
225 : ! get the initial guess for the SCF
226 982 : CALL ls_scf_initial_guess(qs_env, ls_scf_env, nonscf)
227 :
228 982 : 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 982 : 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 982 : 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 982 : CALL timestop(handle)
241 :
242 982 : 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 1452 : 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 470 : CALL timeset(routineN, handle)
273 :
274 : ! get a useful output_unit
275 470 : logger => cp_get_default_logger()
276 470 : IF (logger%para_env%is_source()) THEN
277 235 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
278 : ELSE
279 : unit_nr = -1
280 : END IF
281 :
282 235 : IF (unit_nr > 0) WRITE (unit_nr, '()')
283 : ! if there is no history go for the atomic guess, otherwise extrapolate the dm history
284 470 : IF (ls_scf_env%scf_history%istore == 0) THEN
285 294 : 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 290 : CALL ls_scf_qs_atomic_guess(qs_env, ls_scf_env, ls_scf_env%energy_init, nonscf)
298 290 : 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 176 : CALL cite_reference(Kolafa2004)
319 176 : CALL cite_reference(Kuhne2007)
320 176 : naspc = MIN(ls_scf_env%scf_history%istore, ls_scf_env%scf_history%nstore)
321 358 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
322 : ! actual extrapolation
323 182 : CALL dbcsr_set(ls_scf_env%matrix_p(ispin), 0.0_dp)
324 900 : DO iaspc = 1, naspc
325 : alpha = (-1.0_dp)**(iaspc + 1)*REAL(iaspc, KIND=dp)* &
326 542 : binomial(2*naspc, naspc - iaspc)/binomial(2*naspc - 2, naspc - 1)
327 542 : istore = MOD(ls_scf_env%scf_history%istore - iaspc, ls_scf_env%scf_history%nstore) + 1
328 724 : 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 176 : 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 176 : IF (.NOT. (ls_scf_env%do_pexsi)) THEN
341 358 : 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 182 : 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 182 : CALL dbcsr_filter(ls_scf_env%matrix_p(ispin), ls_scf_env%eps_filter**(2.0_dp/3.0_dp))
346 182 : CALL purify_mcweeny(ls_scf_env%matrix_p(ispin:ispin), ls_scf_env%eps_filter, 3)
347 182 : IF (SIZE(ls_scf_env%matrix_p) == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
348 :
349 358 : 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 182 : 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 182 : 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 182 : filter_eps=ls_scf_env%eps_filter)
358 182 : CALL dbcsr_release(matrix_tmp1)
359 :
360 182 : IF (ls_scf_env%has_s_preconditioner) THEN
361 : CALL apply_matrix_preconditioner(ls_scf_env%matrix_p(ispin), "forward", &
362 176 : 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 646 : IF (nonscf) THEN
370 60 : CALL ls_nonscf_ks(qs_env, ls_scf_env, ls_scf_env%energy_init)
371 : ELSE
372 116 : 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 470 : IF (unit_nr > 0) THEN
377 235 : WRITE (unit_nr, '(T2,A,F20.9)') "Energy with the initial guess:", ls_scf_env%energy_init
378 235 : WRITE (unit_nr, '()')
379 : END IF
380 :
381 470 : CALL timestop(handle)
382 :
383 982 : 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 470 : 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 470 : CALL timeset(routineN, handle)
404 :
405 : ! get a useful output_unit
406 470 : logger => cp_get_default_logger()
407 470 : IF (logger%para_env%is_source()) THEN
408 235 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
409 : ELSE
410 : unit_nr = -1
411 : END IF
412 :
413 470 : 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 470 : IF (ls_scf_env%scf_history%nstore > 0) THEN
432 462 : ls_scf_env%scf_history%istore = ls_scf_env%scf_history%istore + 1
433 952 : DO ispin = 1, SIZE(ls_scf_env%matrix_p)
434 490 : istore = MOD(ls_scf_env%scf_history%istore - 1, ls_scf_env%scf_history%nstore) + 1
435 490 : 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 952 : 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 486 : matrix_type=dbcsr_type_no_symmetry)
443 :
444 486 : IF (ls_scf_env%has_s_preconditioner) THEN
445 : CALL apply_matrix_preconditioner(ls_scf_env%scf_history%matrix(ispin, istore), "backward", &
446 440 : 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 486 : 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 486 : filter_eps=ls_scf_env%eps_filter)
453 486 : CALL dbcsr_release(matrix_tmp1)
454 : END IF
455 :
456 : END DO
457 : END IF
458 :
459 470 : CALL timestop(handle)
460 :
461 470 : 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 982 : 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 982 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
488 982 : 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 982 : CALL timeset(routineN, handle)
494 :
495 : ! get a useful output_unit
496 982 : logger => cp_get_default_logger()
497 982 : IF (logger%para_env%is_source()) THEN
498 491 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
499 : ELSE
500 491 : unit_nr = -1
501 : END IF
502 :
503 982 : nspin = ls_scf_env%nspins
504 :
505 : ! old quantities, useful for mixing
506 5948 : ALLOCATE (matrix_mixing_old(nspin), matrix_ks_deviation(nspin))
507 1992 : DO ispin = 1, nspin
508 1010 : CALL dbcsr_create(matrix_mixing_old(ispin), template=ls_scf_env%matrix_ks(ispin))
509 :
510 1010 : CALL dbcsr_create(matrix_ks_deviation(ispin), template=ls_scf_env%matrix_ks(ispin))
511 1992 : CALL dbcsr_set(matrix_ks_deviation(ispin), 0.0_dp)
512 : END DO
513 2946 : ls_scf_env%homo_spin(:) = 0.0_dp
514 2946 : ls_scf_env%lumo_spin(:) = 0.0_dp
515 :
516 982 : transm_scf_converged = .FALSE.
517 982 : transm_maxscf_reached = .FALSE.
518 :
519 982 : energy_old = 0.0_dp
520 982 : IF (ls_scf_env%scf_history%istore > 0) energy_old = ls_scf_env%energy_init
521 982 : check_convergence = .TRUE.
522 982 : iscf = 0
523 982 : 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 982 : CALL get_qs_env(qs_env, transport_env=transport_env, do_transport=do_transport)
536 :
537 : ! the real SCF loop
538 3156 : DO
539 :
540 : ! check on max SCF or timing/exit
541 3156 : CALL external_control(should_stop, "SCF", start_time=qs_env%start_time, target_time=qs_env%target_time)
542 3156 : 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 == 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 3156 : 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 3110 : t1 = m_walltime()
570 3110 : iscf = iscf + 1
571 :
572 : ! first get a copy of the current KS matrix
573 3110 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
574 6356 : DO ispin = 1, nspin
575 : CALL matrix_qs_to_ls(ls_scf_env%matrix_ks(ispin), matrix_ks(ispin)%matrix, &
576 3246 : ls_scf_env%ls_mstruct, covariant=.TRUE.)
577 3246 : IF (ls_scf_env%has_s_preconditioner) THEN
578 : CALL apply_matrix_preconditioner(ls_scf_env%matrix_ks(ispin), "forward", &
579 1756 : ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
580 : END IF
581 6356 : 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 3110 : 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 6164 : DO ispin = 1, nspin
589 3144 : IF (nonscf) THEN
590 66 : CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin))
591 3078 : ELSE IF (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 3078 : IF (iscf == 1) THEN
595 : ! initialize the mixing matrix with the current state if needed
596 934 : CALL dbcsr_copy(matrix_mixing_old(ispin), ls_scf_env%matrix_ks(ispin))
597 : ELSE
598 2144 : IF (ls_scf_env%ls_diis) THEN ! ------- IF-DIIS+MIX--- START
599 28 : IF (diis_step .AND. (iscf - 1) >= 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 2116 : IF (unit_nr > 0) THEN
634 : WRITE (unit_nr, '(A57)') &
635 1058 : "*********************************************************"
636 : WRITE (unit_nr, '(A23,F5.3,A25,I3)') &
637 1058 : " Using MIXING_FRACTION=", ls_scf_env%mixing_fraction, &
638 2116 : " to mix KS matrix: iscf=", iscf
639 : WRITE (unit_nr, '(A7,F5.3,A6,F5.3,A7)') &
640 1058 : " KS_nw=", ls_scf_env%mixing_fraction, "*KS + ", &
641 2116 : 1.0_dp - ls_scf_env%mixing_fraction, "*KS_old"
642 : WRITE (unit_nr, '(A57)') &
643 1058 : "*********************************************************"
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 2116 : 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 3144 : nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
657 3144 : IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
658 :
659 3144 : IF (do_transport) THEN
660 0 : IF (ls_scf_env%has_s_preconditioner) THEN
661 0 : CPABORT("NOT YET IMPLEMENTED with S preconditioner. ")
662 : END IF
663 0 : IF (ls_scf_env%ls_mstruct%cluster_type /= ls_cluster_atomic) THEN
664 0 : CPABORT("NOT YET IMPLEMENTED with molecular clustering. ")
665 : END IF
666 :
667 0 : extra_scf = maxscf_reached .OR. scf_converged
668 : ! get the current Kohn-Sham matrix (ks) and return matrix_p evaluated using an external C routine
669 : CALL external_scf_method(transport_env, ls_scf_env%matrix_s, matrix_mixing_old(ispin), &
670 : ls_scf_env%matrix_p(ispin), nelectron_spin_real, ls_scf_env%natoms, &
671 0 : energy_diff, iscf, extra_scf)
672 :
673 : ELSE
674 4186 : SELECT CASE (ls_scf_env%purification_method)
675 : CASE (ls_scf_sign)
676 : CALL density_matrix_sign(ls_scf_env%matrix_p(ispin), ls_scf_env%mu_spin(ispin), ls_scf_env%fixed_mu, &
677 : ls_scf_env%sign_method, ls_scf_env%sign_order, matrix_mixing_old(ispin), &
678 : ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, nelectron_spin_real, &
679 : ls_scf_env%eps_filter, ls_scf_env%sign_symmetric, ls_scf_env%submatrix_sign_method, &
680 1042 : ls_scf_env%matrix_s_sqrt_inv)
681 : CASE (ls_scf_tc2)
682 : CALL density_matrix_tc2(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, &
683 : nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), &
684 : ls_scf_env%lumo_spin(ispin), non_monotonic=ls_scf_env%non_monotonic, &
685 : eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
686 284 : iounit=-1)
687 : CASE (ls_scf_trs4)
688 : CALL density_matrix_trs4(ls_scf_env%matrix_p(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s_sqrt_inv, &
689 : nelectron_spin_real, ls_scf_env%eps_filter, ls_scf_env%homo_spin(ispin), &
690 : ls_scf_env%lumo_spin(ispin), ls_scf_env%mu_spin(ispin), &
691 : dynamic_threshold=ls_scf_env%dynamic_threshold, &
692 : matrix_ks_deviation=matrix_ks_deviation(ispin), &
693 : eps_lanczos=ls_scf_env%eps_lanczos, max_iter_lanczos=ls_scf_env%max_iter_lanczos, &
694 1818 : iounit=-1)
695 : CASE (ls_scf_pexsi)
696 0 : IF (ls_scf_env%has_s_preconditioner) THEN
697 0 : CPABORT("S preconditioning not implemented in combination with the PEXSI library. ")
698 : END IF
699 0 : IF (ls_scf_env%ls_mstruct%cluster_type /= ls_cluster_atomic) THEN
700 : CALL cp_abort(__LOCATION__, &
701 0 : "Molecular clustering not implemented in combination with the PEXSI library. ")
702 : END IF
703 : CALL density_matrix_pexsi(ls_scf_env%pexsi, ls_scf_env%matrix_p(ispin), ls_scf_env%pexsi%matrix_w(ispin), &
704 : ls_scf_env%pexsi%kTS(ispin), matrix_mixing_old(ispin), ls_scf_env%matrix_s, &
705 3144 : nelectron_spin_real, ls_scf_env%mu_spin(ispin), iscf, ispin)
706 : END SELECT
707 : END IF
708 :
709 3144 : IF (ls_scf_env%has_s_preconditioner) THEN
710 : CALL apply_matrix_preconditioner(ls_scf_env%matrix_p(ispin), "forward", &
711 1756 : ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
712 : END IF
713 3144 : CALL dbcsr_filter(ls_scf_env%matrix_p(ispin), ls_scf_env%eps_filter)
714 :
715 6164 : IF (ls_scf_env%nspins == 1) CALL dbcsr_scale(ls_scf_env%matrix_p(ispin), 2.0_dp)
716 :
717 : END DO
718 : END IF
719 :
720 : ! compute the corresponding new energy KS matrix and new energy
721 3110 : IF (nonscf) THEN
722 66 : CALL ls_nonscf_energy(qs_env, ls_scf_env)
723 : ELSE
724 3044 : CALL ls_scf_dm_to_ks(qs_env, ls_scf_env, energy_new, iscf)
725 : END IF
726 :
727 3110 : IF (ls_scf_env%do_pexsi) THEN
728 0 : CALL pexsi_to_qs(ls_scf_env, qs_env, kTS=ls_scf_env%pexsi%kTS)
729 : END IF
730 :
731 3110 : t2 = m_walltime()
732 3110 : IF (nonscf) THEN
733 66 : tdiag = t2 - t1
734 66 : CALL qs_nonscf_print_summary(qs_env, tdiag, ls_scf_env%nelectron_total, unit_nr)
735 66 : EXIT
736 : ELSE
737 : ! report current SCF loop
738 3044 : energy_diff = energy_new - energy_old
739 3044 : energy_old = energy_new
740 3044 : IF (unit_nr > 0) THEN
741 1522 : WRITE (unit_nr, *)
742 1522 : WRITE (unit_nr, '(T2,A,I6,F20.9,F20.9,F12.6)') "SCF", iscf, energy_new, energy_diff, t2 - t1
743 1522 : WRITE (unit_nr, *)
744 1522 : CALL m_flush(unit_nr)
745 : END IF
746 : END IF
747 :
748 3044 : IF (do_transport) THEN
749 0 : scf_converged = check_convergence .AND. ABS(energy_diff) < ls_scf_env%eps_scf*ls_scf_env%nelectron_total
750 : ! one extra scf step for post-processing in transmission calculations
751 0 : IF (transport_env%params%method == transport_transmission) THEN
752 0 : IF (transm_scf_converged) EXIT
753 : transm_scf_converged = scf_converged
754 : ELSE
755 0 : IF (scf_converged) THEN
756 0 : IF (unit_nr > 0) WRITE (unit_nr, '(/,T2,A,I5,A/)') "SCF run converged in ", iscf, " steps."
757 : EXIT
758 : END IF
759 : END IF
760 : ELSE
761 : ! exit criterion on the energy only for the time being
762 3044 : IF (check_convergence .AND. ABS(energy_diff) < ls_scf_env%eps_scf*ls_scf_env%nelectron_total) THEN
763 870 : IF (unit_nr > 0) WRITE (unit_nr, '(/,T2,A,I5,A/)') "SCF run converged in ", iscf, " steps."
764 : ! Skip Harris functional calculation if ground-state is NOT converged
765 870 : IF (qs_env%energy_correction) THEN
766 20 : CALL get_qs_env(qs_env, ec_env=ec_env)
767 20 : IF (ec_env%skip_ec) ec_env%do_skip = .FALSE.
768 : END IF
769 : EXIT
770 : END IF
771 : END IF
772 :
773 2174 : IF (ls_scf_env%ls_diis) THEN
774 : ! diis_buffer, buffer with 1) Kohn-Sham history matrix,
775 : ! 2) KS error history matrix (f=KPS-SPK),
776 : ! 3) B matrix (for finding DIIS weighting coefficients)
777 : CALL qs_diis_b_step_4lscf(diis_buffer, qs_env, ls_scf_env, unit_nr, &
778 : iscf, diis_step, eps_diis, nmixing, matrix_s(1)%matrix, &
779 18 : ls_scf_env%eps_filter)
780 : END IF
781 :
782 2174 : IF (ls_scf_env%do_pexsi) THEN
783 : CALL pexsi_set_convergence_tolerance(ls_scf_env%pexsi, energy_diff, &
784 : ls_scf_env%eps_scf*ls_scf_env%nelectron_total, &
785 : ! initialize in second scf step of first SCF cycle:
786 : (iscf == 2) .AND. (ls_scf_env%scf_history%istore == 0), &
787 0 : check_convergence)
788 : END IF
789 :
790 : END DO
791 :
792 : ! free storage
793 982 : IF (ls_scf_env%ls_diis) THEN
794 4 : CALL qs_diis_b_release_sparse(diis_buffer)
795 4 : DEALLOCATE (diis_buffer)
796 : END IF
797 1992 : DO ispin = 1, nspin
798 1010 : CALL dbcsr_release(matrix_mixing_old(ispin))
799 1992 : CALL dbcsr_release(matrix_ks_deviation(ispin))
800 : END DO
801 982 : DEALLOCATE (matrix_mixing_old, matrix_ks_deviation)
802 :
803 982 : CALL timestop(handle)
804 :
805 982 : END SUBROUTINE ls_scf_main
806 :
807 : ! **************************************************************************************************
808 : !> \brief after SCF we have a density matrix, and the self consistent KS matrix
809 : !> analyze its properties.
810 : !> \param qs_env ...
811 : !> \param ls_scf_env ...
812 : !> \par History
813 : !> 2010.10 created [Joost VandeVondele]
814 : !> \author Joost VandeVondele
815 : ! **************************************************************************************************
816 982 : SUBROUTINE ls_scf_post(qs_env, ls_scf_env)
817 : TYPE(qs_environment_type), POINTER :: qs_env
818 : TYPE(ls_scf_env_type) :: ls_scf_env
819 :
820 : CHARACTER(len=*), PARAMETER :: routineN = 'ls_scf_post'
821 :
822 : INTEGER :: handle, ispin, unit_nr
823 : REAL(KIND=dp) :: occ
824 : TYPE(cp_logger_type), POINTER :: logger
825 982 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
826 : TYPE(dft_control_type), POINTER :: dft_control
827 :
828 982 : CALL timeset(routineN, handle)
829 :
830 982 : CALL get_qs_env(qs_env, dft_control=dft_control)
831 :
832 : ! get a useful output_unit
833 982 : logger => cp_get_default_logger()
834 982 : IF (logger%para_env%is_source()) THEN
835 491 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
836 : ELSE
837 491 : unit_nr = -1
838 : END IF
839 :
840 : ! store the matrix for a next scf run
841 982 : IF (.NOT. ls_scf_env%do_pao) THEN
842 470 : CALL ls_scf_store_result(ls_scf_env)
843 : END IF
844 :
845 : ! write homo and lumo energy and occupation (if not already part of the output)
846 982 : IF (ls_scf_env%curvy_steps) THEN
847 18 : CALL post_scf_homo_lumo(ls_scf_env)
848 :
849 : ! always report P occ
850 18 : IF (unit_nr > 0) WRITE (unit_nr, *) ""
851 38 : DO ispin = 1, ls_scf_env%nspins
852 20 : occ = dbcsr_get_occupation(ls_scf_env%matrix_p(ispin))
853 38 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Density matrix (P) occupation ", occ
854 : END DO
855 : END IF
856 :
857 : ! compute the matrix_w if associated
858 982 : IF (ls_scf_env%calculate_forces) THEN
859 194 : CALL get_qs_env(qs_env, matrix_w=matrix_w)
860 194 : CPASSERT(ASSOCIATED(matrix_w))
861 194 : IF (ls_scf_env%do_pexsi) THEN
862 0 : CALL pexsi_to_qs(ls_scf_env, qs_env, matrix_w=ls_scf_env%pexsi%matrix_w)
863 : ELSE
864 194 : CALL calculate_w_matrix_ls(matrix_w, ls_scf_env)
865 : END IF
866 : END IF
867 :
868 : ! compute properties
869 :
870 982 : IF (ls_scf_env%perform_mu_scan) CALL post_scf_mu_scan(ls_scf_env)
871 :
872 982 : IF (ls_scf_env%report_all_sparsities) CALL post_scf_sparsities(ls_scf_env)
873 :
874 982 : IF (dft_control%qs_control%dftb) THEN
875 54 : CALL scf_post_calculation_tb(qs_env, "DFTB", .TRUE.)
876 928 : ELSE IF (dft_control%qs_control%xtb) THEN
877 94 : CALL scf_post_calculation_tb(qs_env, "xTB", .TRUE.)
878 : ELSE
879 834 : CALL write_mo_free_results(qs_env)
880 : END IF
881 :
882 982 : IF (ls_scf_env%chebyshev%compute_chebyshev) CALL compute_chebyshev(qs_env, ls_scf_env)
883 :
884 982 : IF (.TRUE.) CALL post_scf_experiment()
885 :
886 982 : IF (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb) THEN
887 : !
888 : ELSE
889 834 : CALL qs_scf_post_moments(qs_env%input, logger, qs_env, unit_nr)
890 : END IF
891 :
892 : ! clean up used data
893 :
894 982 : CALL dbcsr_release(ls_scf_env%matrix_s)
895 982 : CALL deallocate_curvy_data(ls_scf_env%curvy_data)
896 :
897 982 : IF (ls_scf_env%has_s_preconditioner) THEN
898 426 : CALL dbcsr_release(ls_scf_env%matrix_bs_sqrt)
899 426 : CALL dbcsr_release(ls_scf_env%matrix_bs_sqrt_inv)
900 : END IF
901 :
902 982 : IF (ls_scf_env%needs_s_inv) THEN
903 980 : CALL dbcsr_release(ls_scf_env%matrix_s_inv)
904 : END IF
905 :
906 982 : IF (ls_scf_env%use_s_sqrt) THEN
907 978 : CALL dbcsr_release(ls_scf_env%matrix_s_sqrt)
908 978 : CALL dbcsr_release(ls_scf_env%matrix_s_sqrt_inv)
909 : END IF
910 :
911 1992 : DO ispin = 1, SIZE(ls_scf_env%matrix_ks)
912 1992 : CALL dbcsr_release(ls_scf_env%matrix_ks(ispin))
913 : END DO
914 982 : DEALLOCATE (ls_scf_env%matrix_ks)
915 :
916 982 : IF (ls_scf_env%do_pexsi) THEN
917 0 : CALL pexsi_finalize_scf(ls_scf_env%pexsi, ls_scf_env%mu_spin)
918 : END IF
919 :
920 982 : CALL timestop(handle)
921 :
922 982 : END SUBROUTINE ls_scf_post
923 :
924 : ! **************************************************************************************************
925 : !> \brief Compute the HOMO LUMO energies post SCF
926 : !> \param ls_scf_env ...
927 : !> \par History
928 : !> 2013.06 created [Joost VandeVondele]
929 : !> \author Joost VandeVondele
930 : ! **************************************************************************************************
931 18 : SUBROUTINE post_scf_homo_lumo(ls_scf_env)
932 : TYPE(ls_scf_env_type) :: ls_scf_env
933 :
934 : CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_homo_lumo'
935 :
936 : INTEGER :: handle, ispin, nspin, unit_nr
937 : LOGICAL :: converged
938 : REAL(KIND=dp) :: eps_max, eps_min, homo, lumo
939 : TYPE(cp_logger_type), POINTER :: logger
940 : TYPE(dbcsr_type) :: matrix_k, matrix_p, matrix_tmp
941 :
942 18 : CALL timeset(routineN, handle)
943 :
944 : ! get a useful output_unit
945 18 : logger => cp_get_default_logger()
946 18 : IF (logger%para_env%is_source()) THEN
947 9 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
948 : ELSE
949 9 : unit_nr = -1
950 : END IF
951 :
952 18 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A)') ""
953 :
954 : ! TODO: remove these limitations
955 18 : CPASSERT(.NOT. ls_scf_env%has_s_preconditioner)
956 18 : CPASSERT(ls_scf_env%use_s_sqrt)
957 :
958 18 : nspin = ls_scf_env%nspins
959 :
960 18 : CALL dbcsr_create(matrix_p, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
961 :
962 18 : CALL dbcsr_create(matrix_k, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
963 :
964 18 : CALL dbcsr_create(matrix_tmp, template=ls_scf_env%matrix_p(1), matrix_type=dbcsr_type_no_symmetry)
965 :
966 38 : DO ispin = 1, nspin
967 : ! ortho basis ks
968 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_ks(ispin), &
969 20 : 0.0_dp, matrix_tmp, filter_eps=ls_scf_env%eps_filter)
970 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, ls_scf_env%matrix_s_sqrt_inv, &
971 20 : 0.0_dp, matrix_k, filter_eps=ls_scf_env%eps_filter)
972 :
973 : ! extremal eigenvalues ks
974 : CALL arnoldi_extremal(matrix_k, eps_max, eps_min, max_iter=ls_scf_env%max_iter_lanczos, &
975 20 : threshold=ls_scf_env%eps_lanczos, converged=converged)
976 :
977 : ! ortho basis p
978 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_p(ispin), &
979 20 : 0.0_dp, matrix_tmp, filter_eps=ls_scf_env%eps_filter)
980 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, ls_scf_env%matrix_s_sqrt, &
981 20 : 0.0_dp, matrix_p, filter_eps=ls_scf_env%eps_filter)
982 20 : IF (nspin == 1) CALL dbcsr_scale(matrix_p, 0.5_dp)
983 :
984 : ! go compute homo lumo
985 : CALL compute_homo_lumo(matrix_k, matrix_p, eps_min, eps_max, ls_scf_env%eps_filter, &
986 58 : ls_scf_env%max_iter_lanczos, ls_scf_env%eps_lanczos, homo, lumo, unit_nr)
987 :
988 : END DO
989 :
990 18 : CALL dbcsr_release(matrix_p)
991 18 : CALL dbcsr_release(matrix_k)
992 18 : CALL dbcsr_release(matrix_tmp)
993 :
994 18 : CALL timestop(handle)
995 :
996 18 : END SUBROUTINE post_scf_homo_lumo
997 :
998 : ! **************************************************************************************************
999 : !> \brief Compute the density matrix for various values of the chemical potential
1000 : !> \param ls_scf_env ...
1001 : !> \par History
1002 : !> 2010.10 created [Joost VandeVondele]
1003 : !> \author Joost VandeVondele
1004 : ! **************************************************************************************************
1005 2 : SUBROUTINE post_scf_mu_scan(ls_scf_env)
1006 : TYPE(ls_scf_env_type) :: ls_scf_env
1007 :
1008 : CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_mu_scan'
1009 :
1010 : INTEGER :: handle, imu, ispin, nelectron_spin_real, &
1011 : nmu, nspin, unit_nr
1012 : REAL(KIND=dp) :: mu, t1, t2, trace
1013 : TYPE(cp_logger_type), POINTER :: logger
1014 : TYPE(dbcsr_type) :: matrix_p
1015 :
1016 2 : CALL timeset(routineN, handle)
1017 :
1018 : ! get a useful output_unit
1019 2 : logger => cp_get_default_logger()
1020 2 : IF (logger%para_env%is_source()) THEN
1021 1 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1022 : ELSE
1023 : unit_nr = -1
1024 : END IF
1025 :
1026 2 : nspin = ls_scf_env%nspins
1027 :
1028 2 : CALL dbcsr_create(matrix_p, template=ls_scf_env%matrix_p(1))
1029 :
1030 2 : nmu = 10
1031 24 : DO imu = 0, nmu
1032 :
1033 22 : t1 = m_walltime()
1034 :
1035 22 : mu = -0.4_dp + imu*(0.1_dp + 0.4_dp)/nmu
1036 :
1037 22 : IF (unit_nr > 0) WRITE (unit_nr, *) "------- starting with mu ", mu
1038 :
1039 44 : DO ispin = 1, nspin
1040 : ! we need the proper number of states
1041 22 : nelectron_spin_real = ls_scf_env%nelectron_spin(ispin)
1042 22 : IF (ls_scf_env%nspins == 1) nelectron_spin_real = nelectron_spin_real/2
1043 :
1044 : CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, ls_scf_env%sign_method, &
1045 : ls_scf_env%sign_order, ls_scf_env%matrix_ks(ispin), &
1046 : ls_scf_env%matrix_s, ls_scf_env%matrix_s_inv, &
1047 : ls_scf_env%eps_filter, ls_scf_env%sign_symmetric, &
1048 44 : ls_scf_env%submatrix_sign_method, ls_scf_env%matrix_s_sqrt_inv)
1049 : END DO
1050 :
1051 22 : t2 = m_walltime()
1052 :
1053 24 : IF (unit_nr > 0) WRITE (unit_nr, *) " obtained ", mu, trace, t2 - t1
1054 :
1055 : END DO
1056 :
1057 2 : CALL dbcsr_release(matrix_p)
1058 :
1059 2 : CALL timestop(handle)
1060 :
1061 2 : END SUBROUTINE post_scf_mu_scan
1062 :
1063 : ! **************************************************************************************************
1064 : !> \brief Report on the sparsities of various interesting matrices.
1065 : !>
1066 : !> \param ls_scf_env ...
1067 : !> \par History
1068 : !> 2010.10 created [Joost VandeVondele]
1069 : !> \author Joost VandeVondele
1070 : ! **************************************************************************************************
1071 262 : SUBROUTINE post_scf_sparsities(ls_scf_env)
1072 : TYPE(ls_scf_env_type) :: ls_scf_env
1073 :
1074 : CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_sparsities'
1075 :
1076 : CHARACTER(LEN=default_string_length) :: title
1077 : INTEGER :: handle, ispin, nspin, unit_nr
1078 : TYPE(cp_logger_type), POINTER :: logger
1079 : TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2
1080 :
1081 262 : CALL timeset(routineN, handle)
1082 :
1083 : ! get a useful output_unit
1084 262 : logger => cp_get_default_logger()
1085 262 : IF (logger%para_env%is_source()) THEN
1086 131 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1087 : ELSE
1088 131 : unit_nr = -1
1089 : END IF
1090 :
1091 262 : nspin = ls_scf_env%nspins
1092 :
1093 262 : IF (unit_nr > 0) THEN
1094 131 : WRITE (unit_nr, '()')
1095 131 : WRITE (unit_nr, '(T2,A,E17.3)') "Sparsity reports for eps_filter: ", ls_scf_env%eps_filter
1096 131 : WRITE (unit_nr, '()')
1097 : END IF
1098 :
1099 : CALL report_matrix_sparsity(ls_scf_env%matrix_s, unit_nr, "overlap matrix (S)", &
1100 262 : ls_scf_env%eps_filter)
1101 :
1102 532 : DO ispin = 1, nspin
1103 270 : WRITE (title, '(A,I3)') "Kohn-Sham matrix (H) for spin ", ispin
1104 : CALL report_matrix_sparsity(ls_scf_env%matrix_ks(ispin), unit_nr, title, &
1105 532 : ls_scf_env%eps_filter)
1106 : END DO
1107 :
1108 262 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1109 262 : CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1110 :
1111 532 : DO ispin = 1, nspin
1112 270 : WRITE (title, '(A,I3)') "Density matrix (P) for spin ", ispin
1113 : CALL report_matrix_sparsity(ls_scf_env%matrix_p(ispin), unit_nr, title, &
1114 270 : ls_scf_env%eps_filter)
1115 :
1116 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s, ls_scf_env%matrix_p(ispin), &
1117 270 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1118 :
1119 270 : WRITE (title, '(A,I3,A)') "S * P(", ispin, ")"
1120 270 : CALL report_matrix_sparsity(matrix_tmp1, unit_nr, title, ls_scf_env%eps_filter)
1121 :
1122 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s, &
1123 270 : 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
1124 270 : WRITE (title, '(A,I3,A)') "S * P(", ispin, ") * S"
1125 532 : CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
1126 : END DO
1127 :
1128 262 : IF (ls_scf_env%needs_s_inv) THEN
1129 : CALL report_matrix_sparsity(ls_scf_env%matrix_s_inv, unit_nr, "inv(S)", &
1130 262 : ls_scf_env%eps_filter)
1131 532 : DO ispin = 1, nspin
1132 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_inv, ls_scf_env%matrix_ks(ispin), &
1133 270 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1134 :
1135 270 : WRITE (title, '(A,I3,A)') "inv(S) * H(", ispin, ")"
1136 532 : CALL report_matrix_sparsity(matrix_tmp1, unit_nr, title, ls_scf_env%eps_filter)
1137 : END DO
1138 : END IF
1139 :
1140 262 : IF (ls_scf_env%use_s_sqrt) THEN
1141 :
1142 : CALL report_matrix_sparsity(ls_scf_env%matrix_s_sqrt, unit_nr, "sqrt(S)", &
1143 260 : ls_scf_env%eps_filter)
1144 : CALL report_matrix_sparsity(ls_scf_env%matrix_s_sqrt_inv, unit_nr, "inv(sqrt(S))", &
1145 260 : ls_scf_env%eps_filter)
1146 :
1147 528 : DO ispin = 1, nspin
1148 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_ks(ispin), &
1149 268 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1150 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
1151 268 : 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
1152 268 : WRITE (title, '(A,I3,A)') "inv(sqrt(S)) * H(", ispin, ") * inv(sqrt(S))"
1153 528 : CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
1154 : END DO
1155 :
1156 528 : DO ispin = 1, nspin
1157 : CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_p(ispin), &
1158 268 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1159 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt, &
1160 268 : 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
1161 268 : WRITE (title, '(A,I3,A)') "sqrt(S) * P(", ispin, ") * sqrt(S)"
1162 528 : CALL report_matrix_sparsity(matrix_tmp2, unit_nr, title, ls_scf_env%eps_filter)
1163 : END DO
1164 :
1165 : END IF
1166 :
1167 262 : CALL dbcsr_release(matrix_tmp1)
1168 262 : CALL dbcsr_release(matrix_tmp2)
1169 :
1170 262 : CALL timestop(handle)
1171 :
1172 262 : END SUBROUTINE post_scf_sparsities
1173 :
1174 : ! **************************************************************************************************
1175 : !> \brief Helper routine to report on the sparsity of a single matrix,
1176 : !> for several filtering values
1177 : !> \param matrix ...
1178 : !> \param unit_nr ...
1179 : !> \param title ...
1180 : !> \param eps ...
1181 : !> \par History
1182 : !> 2010.10 created [Joost VandeVondele]
1183 : !> \author Joost VandeVondele
1184 : ! **************************************************************************************************
1185 2930 : SUBROUTINE report_matrix_sparsity(matrix, unit_nr, title, eps)
1186 : TYPE(dbcsr_type) :: matrix
1187 : INTEGER :: unit_nr
1188 : CHARACTER(LEN=*) :: title
1189 : REAL(KIND=dp) :: eps
1190 :
1191 : CHARACTER(len=*), PARAMETER :: routineN = 'report_matrix_sparsity'
1192 :
1193 : INTEGER :: handle
1194 : REAL(KIND=dp) :: eps_local, occ
1195 : TYPE(dbcsr_type) :: matrix_tmp
1196 :
1197 2930 : CALL timeset(routineN, handle)
1198 2930 : CALL dbcsr_create(matrix_tmp, template=matrix, name=TRIM(title))
1199 2930 : CALL dbcsr_copy(matrix_tmp, matrix, name=TRIM(title))
1200 :
1201 2930 : IF (unit_nr > 0) THEN
1202 1465 : WRITE (unit_nr, '(T2,A)') "Sparsity for : "//TRIM(title)
1203 : END IF
1204 :
1205 2930 : eps_local = MAX(eps, 10E-14_dp)
1206 21902 : DO
1207 24832 : IF (eps_local > 1.1_dp) EXIT
1208 21902 : CALL dbcsr_filter(matrix_tmp, eps_local)
1209 21902 : occ = dbcsr_get_occupation(matrix_tmp)
1210 21902 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,F16.12,A3,F16.12)') eps_local, " : ", occ
1211 21902 : eps_local = eps_local*10
1212 : END DO
1213 :
1214 2930 : CALL dbcsr_release(matrix_tmp)
1215 :
1216 2930 : CALL timestop(handle)
1217 :
1218 2930 : END SUBROUTINE report_matrix_sparsity
1219 :
1220 : ! **************************************************************************************************
1221 : !> \brief Compute matrix_w as needed for the forces
1222 : !> \param matrix_w ...
1223 : !> \param ls_scf_env ...
1224 : !> \par History
1225 : !> 2010.11 created [Joost VandeVondele]
1226 : !> \author Joost VandeVondele
1227 : ! **************************************************************************************************
1228 224 : SUBROUTINE calculate_w_matrix_ls(matrix_w, ls_scf_env)
1229 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
1230 : TYPE(ls_scf_env_type) :: ls_scf_env
1231 :
1232 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_ls'
1233 :
1234 : INTEGER :: handle, ispin
1235 : REAL(KIND=dp) :: scaling
1236 : TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2, matrix_tmp3
1237 :
1238 224 : CALL timeset(routineN, handle)
1239 :
1240 224 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1241 224 : CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1242 224 : CALL dbcsr_create(matrix_tmp3, template=ls_scf_env%matrix_s, matrix_type=dbcsr_type_no_symmetry)
1243 :
1244 224 : IF (ls_scf_env%nspins == 1) THEN
1245 216 : scaling = 0.5_dp
1246 : ELSE
1247 8 : scaling = 1.0_dp
1248 : END IF
1249 :
1250 456 : DO ispin = 1, ls_scf_env%nspins
1251 :
1252 232 : CALL dbcsr_copy(matrix_tmp3, ls_scf_env%matrix_ks(ispin))
1253 232 : IF (ls_scf_env%has_s_preconditioner) THEN
1254 : CALL apply_matrix_preconditioner(matrix_tmp3, "backward", &
1255 160 : ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
1256 : END IF
1257 232 : CALL dbcsr_filter(matrix_tmp3, ls_scf_env%eps_filter)
1258 :
1259 : CALL dbcsr_multiply("N", "N", scaling, ls_scf_env%matrix_p(ispin), matrix_tmp3, &
1260 232 : 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
1261 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_p(ispin), &
1262 232 : 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
1263 456 : CALL matrix_ls_to_qs(matrix_w(ispin)%matrix, matrix_tmp2, ls_scf_env%ls_mstruct, covariant=.FALSE.)
1264 : END DO
1265 :
1266 224 : CALL dbcsr_release(matrix_tmp1)
1267 224 : CALL dbcsr_release(matrix_tmp2)
1268 224 : CALL dbcsr_release(matrix_tmp3)
1269 :
1270 224 : CALL timestop(handle)
1271 :
1272 224 : END SUBROUTINE calculate_w_matrix_ls
1273 :
1274 : ! **************************************************************************************************
1275 : !> \brief a place for quick experiments
1276 : !> \par History
1277 : !> 2010.11 created [Joost VandeVondele]
1278 : !> \author Joost VandeVondele
1279 : ! **************************************************************************************************
1280 982 : SUBROUTINE post_scf_experiment()
1281 :
1282 : CHARACTER(len=*), PARAMETER :: routineN = 'post_scf_experiment'
1283 :
1284 : INTEGER :: handle, unit_nr
1285 : TYPE(cp_logger_type), POINTER :: logger
1286 :
1287 982 : CALL timeset(routineN, handle)
1288 :
1289 : ! get a useful output_unit
1290 982 : logger => cp_get_default_logger()
1291 982 : IF (logger%para_env%is_source()) THEN
1292 491 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1293 : ELSE
1294 : unit_nr = -1
1295 : END IF
1296 :
1297 982 : CALL timestop(handle)
1298 982 : END SUBROUTINE post_scf_experiment
1299 :
1300 : END MODULE dm_ls_scf
|
