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