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