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