(git:374b731)
Loading...
Searching...
No Matches
dm_ls_scf_methods.F
Go to the documentation of this file.
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 lower level routines for linear scaling SCF
10!> \par History
11!> 2010.10 created [Joost VandeVondele]
12!> \author Joost VandeVondele
13! **************************************************************************************************
14MODULE dm_ls_scf_methods
15 USE arnoldi_api, ONLY: arnoldi_extremal
16 USE cp_log_handling, ONLY: cp_get_default_logger,&
17 cp_logger_get_default_unit_nr,&
18 cp_logger_type
19 USE dbcsr_api, ONLY: &
20 dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_dot, &
21 dbcsr_filter, dbcsr_finalize, dbcsr_frobenius_norm, dbcsr_get_data_type, &
22 dbcsr_get_occupation, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
23 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, &
24 dbcsr_put_block, dbcsr_release, dbcsr_scale, dbcsr_set, dbcsr_trace, dbcsr_type, &
25 dbcsr_type_no_symmetry
26 USE dm_ls_scf_qs, ONLY: matrix_qs_to_ls
27 USE dm_ls_scf_types, ONLY: ls_cluster_atomic,&
28 ls_mstruct_type,&
29 ls_scf_env_type
30 USE input_constants, ONLY: &
31 ls_cluster_atomic, ls_s_preconditioner_atomic, ls_s_preconditioner_molecular, &
32 ls_s_preconditioner_none, ls_s_sqrt_ns, ls_s_sqrt_proot, ls_scf_sign_ns, &
33 ls_scf_sign_proot, ls_scf_sign_submatrix, ls_scf_submatrix_sign_direct_muadj, &
34 ls_scf_submatrix_sign_direct_muadj_lowmem, ls_scf_submatrix_sign_ns
35 USE iterate_matrix, ONLY: invert_hotelling,&
36 matrix_sign_newton_schulz,&
37 matrix_sign_proot,&
38 matrix_sign_submatrix,&
39 matrix_sign_submatrix_mu_adjust,&
40 matrix_sqrt_newton_schulz,&
41 matrix_sqrt_proot
42 USE kinds, ONLY: dp,&
43 int_8
44 USE machine, ONLY: m_flush,&
45 m_walltime
46 USE mathlib, ONLY: abnormal_value
47#include "./base/base_uses.f90"
48
49 IMPLICIT NONE
50
51 PRIVATE
52
53 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'dm_ls_scf_methods'
54
55 PUBLIC :: ls_scf_init_matrix_s
56 PUBLIC :: density_matrix_sign, density_matrix_sign_fixed_mu
57 PUBLIC :: apply_matrix_preconditioner, compute_matrix_preconditioner
58 PUBLIC :: density_matrix_trs4, density_matrix_tc2, compute_homo_lumo
59
60CONTAINS
61
62! **************************************************************************************************
63!> \brief initialize S matrix related properties (sqrt, inverse...)
64!> Might be factored-out since this seems common code with the other SCF.
65!> \param matrix_s ...
66!> \param ls_scf_env ...
67!> \par History
68!> 2010.10 created [Joost VandeVondele]
69!> \author Joost VandeVondele
70! **************************************************************************************************
71 SUBROUTINE ls_scf_init_matrix_s(matrix_s, ls_scf_env)
72 TYPE(dbcsr_type) :: matrix_s
73 TYPE(ls_scf_env_type) :: ls_scf_env
74
75 CHARACTER(len=*), PARAMETER :: routinen = 'ls_scf_init_matrix_S'
76
77 INTEGER :: handle, unit_nr
78 REAL(kind=dp) :: frob_matrix, frob_matrix_base
79 TYPE(cp_logger_type), POINTER :: logger
80 TYPE(dbcsr_type) :: matrix_tmp1, matrix_tmp2
81
82 CALL timeset(routinen, handle)
83
84 ! get a useful output_unit
85 logger => cp_get_default_logger()
86 IF (logger%para_env%is_source()) THEN
87 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
88 ELSE
89 unit_nr = -1
90 END IF
91
92 ! make our own copy of S
93 IF (ls_scf_env%has_unit_metric) THEN
94 CALL dbcsr_set(ls_scf_env%matrix_s, 0.0_dp)
95 CALL dbcsr_add_on_diag(ls_scf_env%matrix_s, 1.0_dp)
96 ELSE
97 CALL matrix_qs_to_ls(ls_scf_env%matrix_s, matrix_s, ls_scf_env%ls_mstruct, covariant=.true.)
98 END IF
99
100 CALL dbcsr_filter(ls_scf_env%matrix_s, ls_scf_env%eps_filter)
101
102 ! needs a preconditioner for S
103 IF (ls_scf_env%has_s_preconditioner) THEN
104 CALL dbcsr_create(ls_scf_env%matrix_bs_sqrt, template=ls_scf_env%matrix_s, &
105 matrix_type=dbcsr_type_no_symmetry)
106 CALL dbcsr_create(ls_scf_env%matrix_bs_sqrt_inv, template=ls_scf_env%matrix_s, &
107 matrix_type=dbcsr_type_no_symmetry)
108 CALL compute_matrix_preconditioner(ls_scf_env%matrix_s, &
109 ls_scf_env%s_preconditioner_type, ls_scf_env%ls_mstruct, &
110 ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv, &
111 ls_scf_env%eps_filter, ls_scf_env%s_sqrt_order, &
112 ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos)
113 END IF
114
115 ! precondition S
116 IF (ls_scf_env%has_s_preconditioner) THEN
117 CALL apply_matrix_preconditioner(ls_scf_env%matrix_s, "forward", &
118 ls_scf_env%matrix_bs_sqrt, ls_scf_env%matrix_bs_sqrt_inv)
119 END IF
120
121 ! compute sqrt(S) and inv(sqrt(S))
122 IF (ls_scf_env%use_s_sqrt) THEN
123
124 CALL dbcsr_create(ls_scf_env%matrix_s_sqrt, template=ls_scf_env%matrix_s, &
125 matrix_type=dbcsr_type_no_symmetry)
126 CALL dbcsr_create(ls_scf_env%matrix_s_sqrt_inv, template=ls_scf_env%matrix_s, &
127 matrix_type=dbcsr_type_no_symmetry)
128
129 SELECT CASE (ls_scf_env%s_sqrt_method)
130 CASE (ls_s_sqrt_proot)
131 CALL matrix_sqrt_proot(ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_s_sqrt_inv, &
132 ls_scf_env%matrix_s, ls_scf_env%eps_filter, &
133 ls_scf_env%s_sqrt_order, &
134 ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos, &
135 symmetrize=.true.)
136 CASE (ls_s_sqrt_ns)
137 CALL matrix_sqrt_newton_schulz(ls_scf_env%matrix_s_sqrt, ls_scf_env%matrix_s_sqrt_inv, &
138 ls_scf_env%matrix_s, ls_scf_env%eps_filter, &
139 ls_scf_env%s_sqrt_order, &
140 ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos)
141 CASE DEFAULT
142 cpabort("Unknown sqrt method.")
143 END SELECT
144
145 IF (ls_scf_env%check_s_inv) THEN
146 CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
147 matrix_type=dbcsr_type_no_symmetry)
148 CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, &
149 matrix_type=dbcsr_type_no_symmetry)
150
151 CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_s, &
152 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
153
154 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp1, ls_scf_env%matrix_s_sqrt_inv, &
155 0.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
156
157 frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp2)
158 CALL dbcsr_add_on_diag(matrix_tmp2, -1.0_dp)
159 frob_matrix = dbcsr_frobenius_norm(matrix_tmp2)
160 IF (unit_nr > 0) THEN
161 WRITE (unit_nr, *) "Error for (inv(sqrt(S))*S*inv(sqrt(S))-I)", frob_matrix/frob_matrix_base
162 END IF
163
164 CALL dbcsr_release(matrix_tmp1)
165 CALL dbcsr_release(matrix_tmp2)
166 END IF
167 END IF
168
169 ! compute the inverse of S
170 IF (ls_scf_env%needs_s_inv) THEN
171 CALL dbcsr_create(ls_scf_env%matrix_s_inv, template=ls_scf_env%matrix_s, &
172 matrix_type=dbcsr_type_no_symmetry)
173 IF (.NOT. ls_scf_env%use_s_sqrt) THEN
174 CALL invert_hotelling(ls_scf_env%matrix_s_inv, ls_scf_env%matrix_s, ls_scf_env%eps_filter)
175 ELSE
176 CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_sqrt_inv, ls_scf_env%matrix_s_sqrt_inv, &
177 0.0_dp, ls_scf_env%matrix_s_inv, filter_eps=ls_scf_env%eps_filter)
178 END IF
179 IF (ls_scf_env%check_s_inv) THEN
180 CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
181 matrix_type=dbcsr_type_no_symmetry)
182 CALL dbcsr_multiply("N", "N", 1.0_dp, ls_scf_env%matrix_s_inv, ls_scf_env%matrix_s, &
183 0.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
184 frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp1)
185 CALL dbcsr_add_on_diag(matrix_tmp1, -1.0_dp)
186 frob_matrix = dbcsr_frobenius_norm(matrix_tmp1)
187 IF (unit_nr > 0) THEN
188 WRITE (unit_nr, *) "Error for (inv(S)*S-I)", frob_matrix/frob_matrix_base
189 END IF
190 CALL dbcsr_release(matrix_tmp1)
191 END IF
192 END IF
193
194 CALL timestop(handle)
195 END SUBROUTINE ls_scf_init_matrix_s
196
197! **************************************************************************************************
198!> \brief compute for a block positive definite matrix s (bs)
199!> the sqrt(bs) and inv(sqrt(bs))
200!> \param matrix_s ...
201!> \param preconditioner_type ...
202!> \param ls_mstruct ...
203!> \param matrix_bs_sqrt ...
204!> \param matrix_bs_sqrt_inv ...
205!> \param threshold ...
206!> \param order ...
207!> \param eps_lanczos ...
208!> \param max_iter_lanczos ...
209!> \par History
210!> 2010.10 created [Joost VandeVondele]
211!> \author Joost VandeVondele
212! **************************************************************************************************
213 SUBROUTINE compute_matrix_preconditioner(matrix_s, preconditioner_type, ls_mstruct, &
214 matrix_bs_sqrt, matrix_bs_sqrt_inv, threshold, order, eps_lanczos, max_iter_lanczos)
215
216 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_s
217 INTEGER :: preconditioner_type
218 TYPE(ls_mstruct_type) :: ls_mstruct
219 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_bs_sqrt, matrix_bs_sqrt_inv
220 REAL(kind=dp) :: threshold
221 INTEGER :: order
222 REAL(kind=dp) :: eps_lanczos
223 INTEGER :: max_iter_lanczos
224
225 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_matrix_preconditioner'
226
227 INTEGER :: datatype, handle, iblock_col, iblock_row
228 LOGICAL :: block_needed
229 REAL(dp), DIMENSION(:, :), POINTER :: block_dp
230 TYPE(dbcsr_iterator_type) :: iter
231 TYPE(dbcsr_type) :: matrix_bs
232
233 CALL timeset(routinen, handle)
234
235 datatype = dbcsr_get_data_type(matrix_s) ! could be single or double precision
236
237 ! first generate a block diagonal copy of s
238 CALL dbcsr_create(matrix_bs, template=matrix_s)
239
240 SELECT CASE (preconditioner_type)
241 CASE (ls_s_preconditioner_none)
242 CASE (ls_s_preconditioner_atomic, ls_s_preconditioner_molecular)
243 CALL dbcsr_iterator_start(iter, matrix_s)
244 DO WHILE (dbcsr_iterator_blocks_left(iter))
245 CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, block_dp)
246
247 ! do we need the block ?
248 ! this depends on the preconditioner, but also the matrix clustering method employed
249 ! for a clustered matrix, right now, we assume that atomic and molecular preconditioners
250 ! are actually the same, and only require that the diagonal blocks (clustered) are present
251
252 block_needed = .false.
253
254 IF (iblock_row == iblock_col) THEN
255 block_needed = .true.
256 ELSE
257 IF (preconditioner_type == ls_s_preconditioner_molecular .AND. &
258 ls_mstruct%cluster_type == ls_cluster_atomic) THEN
259 IF (ls_mstruct%atom_to_molecule(iblock_row) == ls_mstruct%atom_to_molecule(iblock_col)) block_needed = .true.
260 END IF
261 END IF
262
263 ! add it
264 IF (block_needed) THEN
265 CALL dbcsr_put_block(matrix=matrix_bs, row=iblock_row, col=iblock_col, block=block_dp)
266 END IF
267
268 END DO
269 CALL dbcsr_iterator_stop(iter)
270 END SELECT
271
272 CALL dbcsr_finalize(matrix_bs)
273
274 SELECT CASE (preconditioner_type)
275 CASE (ls_s_preconditioner_none)
276 ! for now make it a simple identity matrix
277 CALL dbcsr_copy(matrix_bs_sqrt, matrix_bs)
278 CALL dbcsr_set(matrix_bs_sqrt, 0.0_dp)
279 CALL dbcsr_add_on_diag(matrix_bs_sqrt, 1.0_dp)
280
281 ! for now make it a simple identity matrix
282 CALL dbcsr_copy(matrix_bs_sqrt_inv, matrix_bs)
283 CALL dbcsr_set(matrix_bs_sqrt_inv, 0.0_dp)
284 CALL dbcsr_add_on_diag(matrix_bs_sqrt_inv, 1.0_dp)
285 CASE (ls_s_preconditioner_atomic, ls_s_preconditioner_molecular)
286 CALL dbcsr_copy(matrix_bs_sqrt, matrix_bs)
287 CALL dbcsr_copy(matrix_bs_sqrt_inv, matrix_bs)
288 ! XXXXXXXXXXX
289 ! XXXXXXXXXXX the threshold here could be done differently,
290 ! XXXXXXXXXXX using eps_filter is reducing accuracy for no good reason, this is cheap
291 ! XXXXXXXXXXX
292 CALL matrix_sqrt_newton_schulz(matrix_bs_sqrt, matrix_bs_sqrt_inv, matrix_bs, &
293 threshold=min(threshold, 1.0e-10_dp), order=order, &
294 eps_lanczos=eps_lanczos, max_iter_lanczos=max_iter_lanczos)
295 END SELECT
296
297 CALL dbcsr_release(matrix_bs)
298
299 CALL timestop(handle)
300
301 END SUBROUTINE compute_matrix_preconditioner
302
303! **************************************************************************************************
304!> \brief apply a preconditioner either
305!> forward (precondition) inv(sqrt(bs)) * A * inv(sqrt(bs))
306!> backward (restore to old form) sqrt(bs) * A * sqrt(bs)
307!> \param matrix ...
308!> \param direction ...
309!> \param matrix_bs_sqrt ...
310!> \param matrix_bs_sqrt_inv ...
311!> \par History
312!> 2010.10 created [Joost VandeVondele]
313!> \author Joost VandeVondele
314! **************************************************************************************************
315 SUBROUTINE apply_matrix_preconditioner(matrix, direction, matrix_bs_sqrt, matrix_bs_sqrt_inv)
316
317 TYPE(dbcsr_type), INTENT(INOUT) :: matrix
318 CHARACTER(LEN=*) :: direction
319 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_bs_sqrt, matrix_bs_sqrt_inv
320
321 CHARACTER(LEN=*), PARAMETER :: routinen = 'apply_matrix_preconditioner'
322
323 INTEGER :: handle
324 TYPE(dbcsr_type) :: matrix_tmp
325
326 CALL timeset(routinen, handle)
327 CALL dbcsr_create(matrix_tmp, template=matrix, matrix_type=dbcsr_type_no_symmetry)
328
329 SELECT CASE (direction)
330 CASE ("forward")
331 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_bs_sqrt_inv, &
332 0.0_dp, matrix_tmp)
333 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_bs_sqrt_inv, matrix_tmp, &
334 0.0_dp, matrix)
335 CASE ("backward")
336 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_bs_sqrt, &
337 0.0_dp, matrix_tmp)
338 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_bs_sqrt, matrix_tmp, &
339 0.0_dp, matrix)
340 CASE DEFAULT
341 cpabort("")
342 END SELECT
343
344 CALL dbcsr_release(matrix_tmp)
345
346 CALL timestop(handle)
347
348 END SUBROUTINE apply_matrix_preconditioner
349
350! **************************************************************************************************
351!> \brief compute the density matrix with a trace that is close to nelectron.
352!> take a mu as input, and improve by bisection as needed.
353!> \param matrix_p ...
354!> \param mu ...
355!> \param fixed_mu ...
356!> \param sign_method ...
357!> \param sign_order ...
358!> \param matrix_ks ...
359!> \param matrix_s ...
360!> \param matrix_s_inv ...
361!> \param nelectron ...
362!> \param threshold ...
363!> \param sign_symmetric ...
364!> \param submatrix_sign_method ...
365!> \param matrix_s_sqrt_inv ...
366!> \par History
367!> 2010.10 created [Joost VandeVondele]
368!> 2020.07 support for methods with internal mu adjustment [Michael Lass]
369!> \author Joost VandeVondele
370! **************************************************************************************************
371 SUBROUTINE density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, &
372 matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric, submatrix_sign_method, &
373 matrix_s_sqrt_inv)
374
375 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
376 REAL(kind=dp), INTENT(INOUT) :: mu
377 LOGICAL :: fixed_mu
378 INTEGER :: sign_method, sign_order
379 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s, matrix_s_inv
380 INTEGER, INTENT(IN) :: nelectron
381 REAL(kind=dp), INTENT(IN) :: threshold
382 LOGICAL, OPTIONAL :: sign_symmetric
383 INTEGER, OPTIONAL :: submatrix_sign_method
384 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_s_sqrt_inv
385
386 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_sign'
387 REAL(kind=dp), PARAMETER :: initial_increment = 0.01_dp
388
389 INTEGER :: handle, iter, unit_nr, &
390 used_submatrix_sign_method
391 LOGICAL :: do_sign_symmetric, has_mu_high, &
392 has_mu_low, internal_mu_adjust
393 REAL(kind=dp) :: increment, mu_high, mu_low, trace
394 TYPE(cp_logger_type), POINTER :: logger
395
396 CALL timeset(routinen, handle)
397
398 logger => cp_get_default_logger()
399 IF (logger%para_env%is_source()) THEN
400 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
401 ELSE
402 unit_nr = -1
403 END IF
404
405 do_sign_symmetric = .false.
406 IF (PRESENT(sign_symmetric)) do_sign_symmetric = sign_symmetric
407
408 used_submatrix_sign_method = ls_scf_submatrix_sign_ns
409 IF (PRESENT(submatrix_sign_method)) used_submatrix_sign_method = submatrix_sign_method
410
411 internal_mu_adjust = ((sign_method .EQ. ls_scf_sign_submatrix) .AND. &
412 (used_submatrix_sign_method .EQ. ls_scf_submatrix_sign_direct_muadj .OR. &
413 used_submatrix_sign_method .EQ. ls_scf_submatrix_sign_direct_muadj_lowmem))
414
415 IF (internal_mu_adjust) THEN
416 CALL density_matrix_sign_internal_mu(matrix_p, trace, mu, sign_method, &
417 matrix_ks, matrix_s, threshold, &
418 used_submatrix_sign_method, &
419 nelectron, matrix_s_sqrt_inv)
420 ELSE
421 increment = initial_increment
422
423 has_mu_low = .false.
424 has_mu_high = .false.
425
426 ! bisect if both bounds are known, otherwise find the bounds with a linear search
427 DO iter = 1, 30
428 IF (has_mu_low .AND. has_mu_high) THEN
429 mu = (mu_low + mu_high)/2
430 IF (abs(mu_high - mu_low) < threshold) EXIT
431 END IF
432
433 CALL density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, &
434 matrix_ks, matrix_s, matrix_s_inv, threshold, &
435 do_sign_symmetric, used_submatrix_sign_method, &
436 matrix_s_sqrt_inv)
437 IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,I2,1X,F13.9,1X,F15.9)') &
438 "Density matrix: iter, mu, trace error: ", iter, mu, trace - nelectron
439
440 ! OK, we can skip early if we are as close as possible to the exact result
441 ! smaller differences should be considered 'noise'
442 IF (abs(trace - nelectron) < 0.5_dp .OR. fixed_mu) EXIT
443
444 IF (trace < nelectron) THEN
445 mu_low = mu
446 mu = mu + increment
447 has_mu_low = .true.
448 increment = increment*2
449 ELSE
450 mu_high = mu
451 mu = mu - increment
452 has_mu_high = .true.
453 increment = increment*2
454 END IF
455 END DO
456
457 END IF
458
459 CALL timestop(handle)
460
461 END SUBROUTINE density_matrix_sign
462
463! **************************************************************************************************
464!> \brief for a fixed mu, compute the corresponding density matrix and its trace
465!> \param matrix_p ...
466!> \param trace ...
467!> \param mu ...
468!> \param sign_method ...
469!> \param sign_order ...
470!> \param matrix_ks ...
471!> \param matrix_s ...
472!> \param matrix_s_inv ...
473!> \param threshold ...
474!> \param sign_symmetric ...
475!> \param submatrix_sign_method ...
476!> \param matrix_s_sqrt_inv ...
477!> \par History
478!> 2010.10 created [Joost VandeVondele]
479!> \author Joost VandeVondele
480! **************************************************************************************************
481 SUBROUTINE density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, matrix_ks, &
482 matrix_s, matrix_s_inv, threshold, sign_symmetric, submatrix_sign_method, &
483 matrix_s_sqrt_inv)
484
485 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
486 REAL(kind=dp), INTENT(OUT) :: trace
487 REAL(kind=dp), INTENT(INOUT) :: mu
488 INTEGER :: sign_method, sign_order
489 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s, matrix_s_inv
490 REAL(kind=dp), INTENT(IN) :: threshold
491 LOGICAL :: sign_symmetric
492 INTEGER :: submatrix_sign_method
493 TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: matrix_s_sqrt_inv
494
495 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_sign_fixed_mu'
496
497 INTEGER :: handle, unit_nr
498 REAL(kind=dp) :: frob_matrix
499 TYPE(cp_logger_type), POINTER :: logger
500 TYPE(dbcsr_type) :: matrix_p_ud, matrix_sign, matrix_sinv_ks, matrix_ssqrtinv_ks_ssqrtinv, &
501 matrix_ssqrtinv_ks_ssqrtinv2, matrix_tmp
502
503 CALL timeset(routinen, handle)
504
505 logger => cp_get_default_logger()
506 IF (logger%para_env%is_source()) THEN
507 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
508 ELSE
509 unit_nr = -1
510 END IF
511
512 CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
513
514 IF (sign_symmetric) THEN
515
516 IF (.NOT. PRESENT(matrix_s_sqrt_inv)) &
517 cpabort("Argument matrix_s_sqrt_inv required if sign_symmetric is set")
518
519 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
520 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv2, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
521 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
522 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, filter_eps=threshold)
523 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, matrix_s_sqrt_inv, &
524 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
525 CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)
526
527 SELECT CASE (sign_method)
528 CASE (ls_scf_sign_ns)
529 CALL matrix_sign_newton_schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
530 CASE (ls_scf_sign_proot)
531 CALL matrix_sign_proot(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
532 CASE (ls_scf_sign_submatrix)
533 CALL matrix_sign_submatrix(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, submatrix_sign_method)
534 CASE DEFAULT
535 cpabort("Unkown sign method.")
536 END SELECT
537 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
538 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2)
539
540 ELSE ! .NOT. sign_symmetric
541 ! get inv(S)*H-I*mu
542 CALL dbcsr_create(matrix_sinv_ks, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
543 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_inv, matrix_ks, &
544 0.0_dp, matrix_sinv_ks, filter_eps=threshold)
545 CALL dbcsr_add_on_diag(matrix_sinv_ks, -mu)
546
547 ! compute sign(inv(S)*H-I*mu)
548 SELECT CASE (sign_method)
549 CASE (ls_scf_sign_ns)
550 CALL matrix_sign_newton_schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order)
551 CASE (ls_scf_sign_proot)
552 CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order)
553 CASE (ls_scf_sign_submatrix)
554 CALL matrix_sign_submatrix(matrix_sign, matrix_sinv_ks, threshold, sign_order, submatrix_sign_method)
555 CASE DEFAULT
556 cpabort("Unkown sign method.")
557 END SELECT
558 CALL dbcsr_release(matrix_sinv_ks)
559 END IF
560
561 ! now construct the density matrix PS=0.5*(I-sign(inv(S)H-I*mu))
562 CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
563 CALL dbcsr_copy(matrix_p_ud, matrix_sign)
564 CALL dbcsr_scale(matrix_p_ud, -0.5_dp)
565 CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp)
566 CALL dbcsr_release(matrix_sign)
567
568 ! we now have PS, lets get its trace
569 CALL dbcsr_trace(matrix_p_ud, trace)
570
571 ! we can also check it is idempotent PS*PS=PS
572 CALL dbcsr_create(matrix_tmp, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
573 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_p_ud, &
574 0.0_dp, matrix_tmp, filter_eps=threshold)
575 CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp)
576 frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
577 IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix
578
579 IF (sign_symmetric) THEN
580 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, &
581 0.0_dp, matrix_tmp, filter_eps=threshold)
582 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_s_sqrt_inv, &
583 0.0_dp, matrix_p, filter_eps=threshold)
584 ELSE
585
586 ! get P=PS*inv(S)
587 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_s_inv, &
588 0.0_dp, matrix_p, filter_eps=threshold)
589 END IF
590 CALL dbcsr_release(matrix_p_ud)
591 CALL dbcsr_release(matrix_tmp)
592
593 CALL timestop(handle)
594
595 END SUBROUTINE density_matrix_sign_fixed_mu
596
597! **************************************************************************************************
598!> \brief compute the corresponding density matrix and its trace, using methods with internal mu adjustment
599!> \param matrix_p ...
600!> \param trace ...
601!> \param mu ...
602!> \param sign_method ...
603!> \param matrix_ks ...
604!> \param matrix_s ...
605!> \param threshold ...
606!> \param submatrix_sign_method ...
607!> \param nelectron ...
608!> \param matrix_s_sqrt_inv ...
609!> \par History
610!> 2020.07 created, based on density_matrix_sign_fixed_mu [Michael Lass]
611!> \author Michael Lass
612! **************************************************************************************************
613 SUBROUTINE density_matrix_sign_internal_mu(matrix_p, trace, mu, sign_method, matrix_ks, &
614 matrix_s, threshold, submatrix_sign_method, &
615 nelectron, matrix_s_sqrt_inv)
616
617 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
618 REAL(kind=dp), INTENT(OUT) :: trace
619 REAL(kind=dp), INTENT(INOUT) :: mu
620 INTEGER :: sign_method
621 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s
622 REAL(kind=dp), INTENT(IN) :: threshold
623 INTEGER :: submatrix_sign_method
624 INTEGER, INTENT(IN) :: nelectron
625 TYPE(dbcsr_type), INTENT(IN) :: matrix_s_sqrt_inv
626
627 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_sign_internal_mu'
628
629 INTEGER :: handle, unit_nr
630 REAL(kind=dp) :: frob_matrix
631 TYPE(cp_logger_type), POINTER :: logger
632 TYPE(dbcsr_type) :: matrix_p_ud, matrix_sign, &
633 matrix_ssqrtinv_ks_ssqrtinv, &
634 matrix_ssqrtinv_ks_ssqrtinv2, &
635 matrix_tmp
636
637 CALL timeset(routinen, handle)
638
639 logger => cp_get_default_logger()
640 IF (logger%para_env%is_source()) THEN
641 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
642 ELSE
643 unit_nr = -1
644 END IF
645
646 CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
647
648 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
649 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv2, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
650 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
651 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, filter_eps=threshold)
652 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, matrix_s_sqrt_inv, &
653 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
654 CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)
655
656 SELECT CASE (sign_method)
657 CASE (ls_scf_sign_submatrix)
658 SELECT CASE (submatrix_sign_method)
659 CASE (ls_scf_submatrix_sign_direct_muadj, ls_scf_submatrix_sign_direct_muadj_lowmem)
660 CALL matrix_sign_submatrix_mu_adjust(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, mu, nelectron, threshold, &
661 submatrix_sign_method)
662 CASE DEFAULT
663 cpabort("density_matrix_sign_internal_mu called with invalid submatrix sign method")
664 END SELECT
665 CASE DEFAULT
666 cpabort("density_matrix_sign_internal_mu called with invalid sign method.")
667 END SELECT
668 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
669 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2)
670
671 ! now construct the density matrix PS=0.5*(I-sign(inv(S)H-I*mu))
672 CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
673 CALL dbcsr_copy(matrix_p_ud, matrix_sign)
674 CALL dbcsr_scale(matrix_p_ud, -0.5_dp)
675 CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp)
676 CALL dbcsr_release(matrix_sign)
677
678 ! we now have PS, lets get its trace
679 CALL dbcsr_trace(matrix_p_ud, trace)
680
681 ! we can also check it is idempotent PS*PS=PS
682 CALL dbcsr_create(matrix_tmp, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
683 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_p_ud, &
684 0.0_dp, matrix_tmp, filter_eps=threshold)
685 CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp)
686 frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
687 IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix
688
689 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, &
690 0.0_dp, matrix_tmp, filter_eps=threshold)
691 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_s_sqrt_inv, &
692 0.0_dp, matrix_p, filter_eps=threshold)
693 CALL dbcsr_release(matrix_p_ud)
694 CALL dbcsr_release(matrix_tmp)
695
696 CALL timestop(handle)
697
698 END SUBROUTINE density_matrix_sign_internal_mu
699
700! **************************************************************************************************
701!> \brief compute the density matrix using a trace-resetting algorithm
702!> \param matrix_p ...
703!> \param matrix_ks ...
704!> \param matrix_s_sqrt_inv ...
705!> \param nelectron ...
706!> \param threshold ...
707!> \param e_homo ...
708!> \param e_lumo ...
709!> \param e_mu ...
710!> \param dynamic_threshold ...
711!> \param matrix_ks_deviation ...
712!> \param max_iter_lanczos ...
713!> \param eps_lanczos ...
714!> \param converged ...
715!> \par History
716!> 2012.06 created [Florian Thoele]
717!> \author Florian Thoele
718! **************************************************************************************************
719 SUBROUTINE density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, &
720 nelectron, threshold, e_homo, e_lumo, e_mu, &
721 dynamic_threshold, matrix_ks_deviation, &
722 max_iter_lanczos, eps_lanczos, converged)
723
724 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
725 TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv
726 INTEGER, INTENT(IN) :: nelectron
727 REAL(kind=dp), INTENT(IN) :: threshold
728 REAL(kind=dp), INTENT(INOUT) :: e_homo, e_lumo, e_mu
729 LOGICAL, INTENT(IN), OPTIONAL :: dynamic_threshold
730 TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: matrix_ks_deviation
731 INTEGER, INTENT(IN) :: max_iter_lanczos
732 REAL(kind=dp), INTENT(IN) :: eps_lanczos
733 LOGICAL, INTENT(OUT), OPTIONAL :: converged
734
735 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_trs4'
736 INTEGER, PARAMETER :: max_iter = 100
737 REAL(kind=dp), PARAMETER :: gamma_max = 6.0_dp, gamma_min = 0.0_dp
738
739 INTEGER :: branch, estimated_steps, handle, i, j, &
740 unit_nr
741 INTEGER(kind=int_8) :: flop1, flop2
742 LOGICAL :: arnoldi_converged, do_dyn_threshold
743 REAL(kind=dp) :: current_threshold, delta_n, eps_max, eps_min, est_threshold, frob_id, &
744 frob_x, gam, homo, lumo, max_eig, max_threshold, maxdev, maxev, min_eig, minev, mmin, mu, &
745 mu_a, mu_b, mu_c, mu_fa, mu_fc, occ_matrix, scaled_homo_bound, scaled_lumo_bound, t1, t2, &
746 trace_fx, trace_gx, xi
747 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: gamma_values
748 TYPE(cp_logger_type), POINTER :: logger
749 TYPE(dbcsr_type) :: matrix_k0, matrix_x, matrix_x_nosym, &
750 matrix_xidsq, matrix_xsq, tmp_gx
751
752 IF (nelectron == 0) THEN
753 CALL dbcsr_set(matrix_p, 0.0_dp)
754 RETURN
755 END IF
756
757 CALL timeset(routinen, handle)
758
759 logger => cp_get_default_logger()
760 IF (logger%para_env%is_source()) THEN
761 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
762 ELSE
763 unit_nr = -1
764 END IF
765
766 do_dyn_threshold = .false.
767 IF (PRESENT(dynamic_threshold)) do_dyn_threshold = dynamic_threshold
768
769 IF (PRESENT(converged)) converged = .false.
770
771 ! init X = (eps_n*I - H)/(eps_n - eps_0) ... H* = S^-1/2*H*S^-1/2
772 CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type="S")
773
774 ! at some points the non-symmetric version of x is required
775 CALL dbcsr_create(matrix_x_nosym, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
776
777 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
778 0.0_dp, matrix_x_nosym, filter_eps=threshold)
779 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x_nosym, matrix_s_sqrt_inv, &
780 0.0_dp, matrix_x, filter_eps=threshold)
781 CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym)
782
783 CALL dbcsr_create(matrix_k0, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
784 CALL dbcsr_copy(matrix_k0, matrix_x_nosym)
785
786 ! compute the deviation in the mixed matrix, as seen in the ortho basis
787 IF (do_dyn_threshold) THEN
788 cpassert(PRESENT(matrix_ks_deviation))
789 CALL dbcsr_add(matrix_ks_deviation, matrix_x_nosym, -1.0_dp, 1.0_dp)
790 CALL arnoldi_extremal(matrix_ks_deviation, maxev, minev, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
791 converged=arnoldi_converged)
792 maxdev = max(abs(maxev), abs(minev))
793 IF (unit_nr > 0) THEN
794 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", arnoldi_converged
795 WRITE (unit_nr, '(T6,A,1X,F12.5)') "change in mixed matrix: ", maxdev
796 WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo + maxdev
797 WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo - maxdev
798 WRITE (unit_nr, '(T6,A,1X,L12)') "Predicts a gap ? ", ((e_lumo - maxdev) - (e_homo + maxdev)) > 0
799 END IF
800 ! save the old mixed matrix
801 CALL dbcsr_copy(matrix_ks_deviation, matrix_x_nosym)
802
803 END IF
804
805 ! get largest/smallest eigenvalues for scaling
806 CALL arnoldi_extremal(matrix_x_nosym, max_eig, min_eig, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
807 converged=arnoldi_converged)
808 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", &
809 min_eig, max_eig, " converged: ", arnoldi_converged
810 eps_max = max_eig
811 eps_min = min_eig
812
813 ! scale KS matrix
814 IF (eps_max == eps_min) THEN
815 CALL dbcsr_scale(matrix_x, 1.0_dp/eps_max)
816 ELSE
817 CALL dbcsr_add_on_diag(matrix_x, -eps_max)
818 CALL dbcsr_scale(matrix_x, -1.0_dp/(eps_max - eps_min))
819 END IF
820
821 current_threshold = threshold
822 IF (do_dyn_threshold) THEN
823 ! scale bounds for HOMO/LUMO
824 scaled_homo_bound = (eps_max - (e_homo + maxdev))/(eps_max - eps_min)
825 scaled_lumo_bound = (eps_max - (e_lumo - maxdev))/(eps_max - eps_min)
826 END IF
827
828 CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type="S")
829
830 CALL dbcsr_create(matrix_xidsq, template=matrix_ks, matrix_type="S")
831
832 CALL dbcsr_create(tmp_gx, template=matrix_ks, matrix_type="S")
833
834 ALLOCATE (gamma_values(max_iter))
835
836 DO i = 1, max_iter
837 t1 = m_walltime()
838 flop1 = 0; flop2 = 0
839
840 ! get X*X
841 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, &
842 0.0_dp, matrix_xsq, &
843 filter_eps=current_threshold, flop=flop1)
844
845 ! intermediate use matrix_xidsq to compute = X*X-X
846 CALL dbcsr_copy(matrix_xidsq, matrix_x)
847 CALL dbcsr_add(matrix_xidsq, matrix_xsq, -1.0_dp, 1.0_dp)
848 frob_id = dbcsr_frobenius_norm(matrix_xidsq)
849 frob_x = dbcsr_frobenius_norm(matrix_x)
850
851 ! xidsq = (1-X)*(1-X)
852 ! use (1-x)*(1-x) = 1 + x*x - 2*x
853 CALL dbcsr_copy(matrix_xidsq, matrix_x)
854 CALL dbcsr_add(matrix_xidsq, matrix_xsq, -2.0_dp, 1.0_dp)
855 CALL dbcsr_add_on_diag(matrix_xidsq, 1.0_dp)
856
857 ! tmp_gx = 4X-3X*X
858 CALL dbcsr_copy(tmp_gx, matrix_x)
859 CALL dbcsr_add(tmp_gx, matrix_xsq, 4.0_dp, -3.0_dp)
860
861 ! get gamma
862 ! Tr(F) = Tr(XX*tmp_gx) Tr(G) is equivalent
863 CALL dbcsr_dot(matrix_xsq, matrix_xidsq, trace_gx)
864 CALL dbcsr_dot(matrix_xsq, tmp_gx, trace_fx)
865
866 ! if converged, and gam becomes noisy, fix it to 3, which results in a final McWeeny step.
867 ! do this only if the electron count is reasonable.
868 ! maybe tune if the current criterion is not good enough
869 delta_n = nelectron - trace_fx
870 ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ...
871 IF (((frob_id*frob_id) < (threshold*frob_x*frob_x)) .AND. (abs(delta_n) < 0.5_dp)) THEN
872 gam = 3.0_dp
873 ELSE IF (abs(delta_n) < 1e-14_dp) THEN
874 gam = 0.0_dp ! rare case of perfect electron count
875 ELSE
876 ! make sure, we don't divide by zero, as soon as gam is outside the interval gam_min,gam_max, it doesn't matter
877 gam = delta_n/max(trace_gx, abs(delta_n)/100)
878 END IF
879 gamma_values(i) = gam
880
881 IF (unit_nr > 0 .AND. .false.) THEN
882 WRITE (unit_nr, *) "trace_fx", trace_fx, "trace_gx", trace_gx, "gam", gam, &
883 "frob_id", frob_id, "conv", abs(frob_id/frob_x)
884 END IF
885
886 IF (do_dyn_threshold) THEN
887 ! quantities used for dynamic thresholding, when the estimated gap is larger than zero
888 xi = (scaled_homo_bound - scaled_lumo_bound)
889 IF (xi > 0.0_dp) THEN
890 mmin = 0.5*(scaled_homo_bound + scaled_lumo_bound)
891 max_threshold = abs(1 - 2*mmin)*xi
892
893 scaled_homo_bound = evaluate_trs4_polynomial(scaled_homo_bound, gamma_values(i:), 1)
894 scaled_lumo_bound = evaluate_trs4_polynomial(scaled_lumo_bound, gamma_values(i:), 1)
895 estimated_steps = estimate_steps(scaled_homo_bound, scaled_lumo_bound, threshold)
896
897 est_threshold = (threshold/(estimated_steps + i + 1))*xi/(1 + threshold/(estimated_steps + i + 1))
898 est_threshold = min(max_threshold, est_threshold)
899 IF (i > 1) est_threshold = max(est_threshold, 0.1_dp*current_threshold)
900 current_threshold = est_threshold
901 ELSE
902 current_threshold = threshold
903 END IF
904 END IF
905
906 IF (gam > gamma_max) THEN
907 ! Xn+1 = 2X-X*X
908 CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp)
909 CALL dbcsr_filter(matrix_x, current_threshold)
910 branch = 1
911 ELSE IF (gam < gamma_min) THEN
912 ! Xn+1 = X*X
913 CALL dbcsr_copy(matrix_x, matrix_xsq)
914 branch = 2
915 ELSE
916 ! Xn+1 = F(X) + gam*G(X)
917 CALL dbcsr_add(tmp_gx, matrix_xidsq, 1.0_dp, gam)
918 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_xsq, tmp_gx, &
919 0.0_dp, matrix_x, &
920 flop=flop2, filter_eps=current_threshold)
921 branch = 3
922 END IF
923
924 occ_matrix = dbcsr_get_occupation(matrix_x)
925 t2 = m_walltime()
926 IF (unit_nr > 0) THEN
927 WRITE (unit_nr, &
928 '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TRS4 it ", &
929 i, occ_matrix, abs(trace_gx), t2 - t1, &
930 (flop1 + flop2)/(1.0e6_dp*max(t2 - t1, 0.001_dp)), current_threshold
931 CALL m_flush(unit_nr)
932 END IF
933
934 IF (abnormal_value(trace_gx)) &
935 cpabort("trace_gx is an abnormal value (NaN/Inf).")
936
937 ! a branch of 1 or 2 appears to lead to a less accurate electron number count and premature exit
938 ! if it turns out this does not exit because we get stuck in branch 1/2 for a reason we need to refine further
939 ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ...
940 IF ((frob_id*frob_id) < (threshold*frob_x*frob_x) .AND. branch == 3 .AND. (abs(delta_n) < 0.5_dp)) THEN
941 IF (PRESENT(converged)) converged = .true.
942 EXIT
943 END IF
944
945 END DO
946
947 occ_matrix = dbcsr_get_occupation(matrix_x)
948 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,I3,1X,F10.8,E12.3)') 'Final TRS4 iteration ', i, occ_matrix, abs(trace_gx)
949
950 ! free some memory
951 CALL dbcsr_release(tmp_gx)
952 CALL dbcsr_release(matrix_xsq)
953 CALL dbcsr_release(matrix_xidsq)
954
955 ! output to matrix_p, P = inv(S)^0.5 X inv(S)^0.5
956 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, &
957 0.0_dp, matrix_x_nosym, filter_eps=threshold)
958 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_x_nosym, &
959 0.0_dp, matrix_p, filter_eps=threshold)
960
961 ! calculate the chemical potential by doing a bisection of fk(x0)-0.5, where fk is evaluated using the stored values for gamma
962 ! E. Rubensson et al., Chem Phys Lett 432, 2006, 591-594
963 mu_a = 0.0_dp; mu_b = 1.0_dp;
964 mu_fa = evaluate_trs4_polynomial(mu_a, gamma_values, i - 1) - 0.5_dp
965 DO j = 1, 40
966 mu_c = 0.5*(mu_a + mu_b)
967 mu_fc = evaluate_trs4_polynomial(mu_c, gamma_values, i - 1) - 0.5_dp ! i-1 because in the last iteration, only convergence is checked
968 IF (abs(mu_fc) < 1.0e-6_dp .OR. (mu_b - mu_a)/2 < 1.0e-6_dp) EXIT !TODO: define threshold values
969
970 IF (mu_fc*mu_fa > 0) THEN
971 mu_a = mu_c
972 mu_fa = mu_fc
973 ELSE
974 mu_b = mu_c
975 END IF
976 END DO
977 mu = (eps_min - eps_max)*mu_c + eps_max
978 DEALLOCATE (gamma_values)
979 IF (unit_nr > 0) THEN
980 WRITE (unit_nr, '(T6,A,1X,F12.5)') 'Chemical potential (mu): ', mu
981 END IF
982 e_mu = mu
983
984 IF (do_dyn_threshold) THEN
985 CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym)
986 CALL compute_homo_lumo(matrix_k0, matrix_x_nosym, eps_min, eps_max, &
987 threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr)
988 e_homo = homo
989 e_lumo = lumo
990 END IF
991
992 CALL dbcsr_release(matrix_x)
993 CALL dbcsr_release(matrix_x_nosym)
994 CALL dbcsr_release(matrix_k0)
995 CALL timestop(handle)
996
997 END SUBROUTINE density_matrix_trs4
998
999! **************************************************************************************************
1000!> \brief compute the density matrix using a non monotonic trace conserving
1001!> algorithm based on SIAM DOI. 10.1137/130911585.
1002!> 2014.04 created [Jonathan Mullin]
1003!> \param matrix_p ...
1004!> \param matrix_ks ...
1005!> \param matrix_s_sqrt_inv ...
1006!> \param nelectron ...
1007!> \param threshold ...
1008!> \param e_homo ...
1009!> \param e_lumo ...
1010!> \param non_monotonic ...
1011!> \param eps_lanczos ...
1012!> \param max_iter_lanczos ...
1013!> \author Jonathan Mullin
1014! **************************************************************************************************
1015 SUBROUTINE density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, &
1016 nelectron, threshold, e_homo, e_lumo, &
1017 non_monotonic, eps_lanczos, max_iter_lanczos)
1018
1019 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
1020 TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv
1021 INTEGER, INTENT(IN) :: nelectron
1022 REAL(kind=dp), INTENT(IN) :: threshold
1023 REAL(kind=dp), INTENT(INOUT) :: e_homo, e_lumo
1024 LOGICAL, INTENT(IN), OPTIONAL :: non_monotonic
1025 REAL(kind=dp), INTENT(IN) :: eps_lanczos
1026 INTEGER, INTENT(IN) :: max_iter_lanczos
1027
1028 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_tc2'
1029 INTEGER, PARAMETER :: max_iter = 100
1030
1031 INTEGER :: handle, i, j, k, unit_nr
1032 INTEGER(kind=int_8) :: flop1, flop2
1033 LOGICAL :: converged, do_non_monotonic
1034 REAL(kind=dp) :: beta, betab, eps_max, eps_min, gama, &
1035 max_eig, min_eig, occ_matrix, t1, t2, &
1036 trace_fx, trace_gx
1037 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: alpha, lambda, nu, poly, wu, x, y
1038 TYPE(cp_logger_type), POINTER :: logger
1039 TYPE(dbcsr_type) :: matrix_tmp, matrix_x, matrix_xsq
1040
1041 CALL timeset(routinen, handle)
1042
1043 logger => cp_get_default_logger()
1044 IF (logger%para_env%is_source()) THEN
1045 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1046 ELSE
1047 unit_nr = -1
1048 END IF
1049
1050 do_non_monotonic = .false.
1051 IF (PRESENT(non_monotonic)) do_non_monotonic = non_monotonic
1052
1053 ! init X = (eps_n*I - H)/(eps_n - eps_0) ... H* = S^-1/2*H*S^-1/2
1054 CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1055 CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1056
1057 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
1058 0.0_dp, matrix_xsq, filter_eps=threshold)
1059 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_xsq, matrix_s_sqrt_inv, &
1060 0.0_dp, matrix_x, filter_eps=threshold)
1061
1062 IF (unit_nr > 0) THEN
1063 WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo
1064 WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo
1065 WRITE (unit_nr, '(T6,A,1X,L12)') "Predicts a gap ? ", ((e_lumo) - (e_homo)) > 0
1066 END IF
1067
1068 ! get largest/smallest eigenvalues for scaling
1069 CALL arnoldi_extremal(matrix_x, max_eig, min_eig, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
1070 converged=converged)
1071 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", &
1072 min_eig, max_eig, " converged: ", converged
1073
1074 eps_max = max_eig
1075 eps_min = min_eig
1076
1077 ! scale KS matrix
1078 CALL dbcsr_scale(matrix_x, -1.0_dp)
1079 CALL dbcsr_add_on_diag(matrix_x, eps_max)
1080 CALL dbcsr_scale(matrix_x, 1/(eps_max - eps_min))
1081
1082 CALL dbcsr_copy(matrix_xsq, matrix_x)
1083
1084 CALL dbcsr_create(matrix_tmp, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1085
1086 ALLOCATE (poly(max_iter))
1087 ALLOCATE (nu(max_iter))
1088 ALLOCATE (wu(max_iter))
1089 ALLOCATE (alpha(max_iter))
1090 ALLOCATE (x(4))
1091 ALLOCATE (y(4))
1092 ALLOCATE (lambda(4))
1093
1094! Controls over the non_monotonic bounds, First if low gap, bias slightly
1095 beta = (eps_max - abs(e_lumo))/(eps_max - eps_min)
1096 betab = (eps_max + abs(e_homo))/(eps_max - eps_min)
1097
1098 IF ((beta - betab) < 0.005_dp) THEN
1099 beta = beta - 0.002_dp
1100 betab = betab + 0.002_dp
1101 END IF
1102! Check if input specifies to use monotonic bounds.
1103 IF (.NOT. do_non_monotonic) THEN
1104 beta = 0.0_dp
1105 betab = 1.0_dp
1106 END IF
1107! initial SCF cycle has no reliable estimate of homo/lumo, force monotinic bounds.
1108 IF (e_homo == 0.0_dp) THEN
1109 beta = 0.0_dp
1110 betab = 1.0_dp
1111 END IF
1112
1113 ! init to take true branch first
1114 trace_fx = nelectron
1115 trace_gx = 0
1116
1117 DO i = 1, max_iter
1118 t1 = m_walltime()
1119 flop1 = 0; flop2 = 0
1120
1121 IF (abs(trace_fx - nelectron) <= abs(trace_gx - nelectron)) THEN
1122! Xn+1 = (aX+ (1-a)I)^2
1123 poly(i) = 1.0_dp
1124 alpha(i) = 2.0_dp/(2.0_dp - beta)
1125
1126 CALL dbcsr_scale(matrix_x, alpha(i))
1127 CALL dbcsr_add_on_diag(matrix_x, 1.0_dp - alpha(i))
1128 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, &
1129 0.0_dp, matrix_xsq, &
1130 filter_eps=threshold, flop=flop1)
1131
1132!save X for control variables
1133 CALL dbcsr_copy(matrix_tmp, matrix_x)
1134
1135 CALL dbcsr_copy(matrix_x, matrix_xsq)
1136
1137 beta = (1.0_dp - alpha(i)) + alpha(i)*beta
1138 beta = beta*beta
1139 betab = (1.0_dp - alpha(i)) + alpha(i)*betab
1140 betab = betab*betab
1141 ELSE
1142! Xn+1 = 2aX-a^2*X^2
1143 poly(i) = 0.0_dp
1144 alpha(i) = 2.0_dp/(1.0_dp + betab)
1145
1146 CALL dbcsr_scale(matrix_x, alpha(i))
1147 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, &
1148 0.0_dp, matrix_xsq, &
1149 filter_eps=threshold, flop=flop1)
1150
1151!save X for control variables
1152 CALL dbcsr_copy(matrix_tmp, matrix_x)
1153!
1154 CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp)
1155
1156 beta = alpha(i)*beta
1157 beta = 2.0_dp*beta - beta*beta
1158 betab = alpha(i)*betab
1159 betab = 2.0_dp*betab - betab*betab
1160
1161 END IF
1162 occ_matrix = dbcsr_get_occupation(matrix_x)
1163 t2 = m_walltime()
1164 IF (unit_nr > 0) THEN
1165 WRITE (unit_nr, &
1166 '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TC2 it ", &
1167 i, occ_matrix, t2 - t1, &
1168 (flop1 + flop2)/(1.0e6_dp*(t2 - t1)), threshold
1169 CALL m_flush(unit_nr)
1170 END IF
1171
1172! calculate control terms
1173 CALL dbcsr_trace(matrix_xsq, trace_fx)
1174
1175! intermediate use matrix_xsq compute X- X*X , temorarily use trace_gx
1176 CALL dbcsr_add(matrix_xsq, matrix_tmp, -1.0_dp, 1.0_dp)
1177 CALL dbcsr_trace(matrix_xsq, trace_gx)
1178 nu(i) = dbcsr_frobenius_norm(matrix_xsq)
1179 wu(i) = trace_gx
1180
1181! intermediate use matrix_xsq to compute = 2X - X*X
1182 CALL dbcsr_add(matrix_xsq, matrix_tmp, 1.0_dp, 1.0_dp)
1183 CALL dbcsr_trace(matrix_xsq, trace_gx)
1184! TC2 has quadratic convergence, using the frobeniums norm as an idempotency deviation test.
1185 IF (abs(nu(i)) < (threshold)) EXIT
1186 END DO
1187
1188 occ_matrix = dbcsr_get_occupation(matrix_x)
1189 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,I3,1X,1F10.8,1X,1F10.8)') 'Final TC2 iteration ', i, occ_matrix, abs(nu(i))
1190
1191 ! output to matrix_p, P = inv(S)^0.5 X inv(S)^0.5
1192 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, &
1193 0.0_dp, matrix_tmp, filter_eps=threshold)
1194 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_tmp, &
1195 0.0_dp, matrix_p, filter_eps=threshold)
1196
1197 CALL dbcsr_release(matrix_xsq)
1198 CALL dbcsr_release(matrix_tmp)
1199
1200 ! ALGO 3 from. SIAM DOI. 10.1137/130911585
1201 x(1) = 1.0_dp
1202 x(2) = 1.0_dp
1203 x(3) = 0.0_dp
1204 x(4) = 0.0_dp
1205 gama = 6.0_dp - 4.0_dp*(sqrt(2.0_dp))
1206 gama = gama - gama*gama
1207 DO WHILE (nu(i) < gama)
1208 ! safeguard against negative root, is skipping correct?
1209 IF (wu(i) < 1.0e-14_dp) THEN
1210 i = i - 1
1211 cycle
1212 END IF
1213 IF ((1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)) < 0.0_dp) THEN
1214 i = i - 1
1215 cycle
1216 END IF
1217 y(1) = 0.5_dp*(1.0_dp - sqrt(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)))
1218 y(2) = 0.5_dp*(1.0_dp - sqrt(1.0_dp - 4.0_dp*nu(i)))
1219 y(3) = 0.5_dp*(1.0_dp + sqrt(1.0_dp - 4.0_dp*nu(i)))
1220 y(4) = 0.5_dp*(1.0_dp + sqrt(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)))
1221 y(:) = min(1.0_dp, max(0.0_dp, y(:)))
1222 DO j = i, 1, -1
1223 IF (poly(j) == 1.0_dp) THEN
1224 DO k = 1, 4
1225 y(k) = sqrt(y(k))
1226 y(k) = (y(k) - 1.0_dp + alpha(j))/alpha(j)
1227 END DO ! end K
1228 ELSE
1229 DO k = 1, 4
1230 y(k) = 1.0_dp - sqrt(1.0_dp - y(k))
1231 y(k) = y(k)/alpha(j)
1232 END DO ! end K
1233 END IF ! end poly
1234 END DO ! end j
1235 x(1) = min(x(1), y(1))
1236 x(2) = min(x(2), y(2))
1237 x(3) = max(x(3), y(3))
1238 x(4) = max(x(4), y(4))
1239 i = i - 1
1240 END DO ! end i
1241! lambda 1,2,3,4 are:: out lumo, in lumo, in homo, out homo
1242 DO k = 1, 4
1243 lambda(k) = eps_max - (eps_max - eps_min)*x(k)
1244 END DO ! end k
1245! END ALGO 3 from. SIAM DOI. 10.1137/130911585
1246 e_homo = lambda(4)
1247 e_lumo = lambda(1)
1248 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "outer homo/lumo/gap", e_homo, e_lumo, (e_lumo - e_homo)
1249 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "inner homo/lumo/gap", lambda(3), lambda(2), (lambda(2) - lambda(3))
1250
1251 DEALLOCATE (poly)
1252 DEALLOCATE (nu)
1253 DEALLOCATE (wu)
1254 DEALLOCATE (alpha)
1255 DEALLOCATE (x)
1256 DEALLOCATE (y)
1257 DEALLOCATE (lambda)
1258
1259 CALL dbcsr_release(matrix_x)
1260 CALL timestop(handle)
1261
1262 END SUBROUTINE density_matrix_tc2
1263
1264! **************************************************************************************************
1265!> \brief compute the homo and lumo given a KS matrix and a density matrix in the orthonormalized basis
1266!> and the eps_min and eps_max, min and max eigenvalue of the ks matrix
1267!> \param matrix_k ...
1268!> \param matrix_p ...
1269!> \param eps_min ...
1270!> \param eps_max ...
1271!> \param threshold ...
1272!> \param max_iter_lanczos ...
1273!> \param eps_lanczos ...
1274!> \param homo ...
1275!> \param lumo ...
1276!> \param unit_nr ...
1277!> \par History
1278!> 2012.06 created [Florian Thoele]
1279!> \author Florian Thoele
1280! **************************************************************************************************
1281 SUBROUTINE compute_homo_lumo(matrix_k, matrix_p, eps_min, eps_max, threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr)
1282 TYPE(dbcsr_type) :: matrix_k, matrix_p
1283 REAL(kind=dp) :: eps_min, eps_max, threshold
1284 INTEGER, INTENT(IN) :: max_iter_lanczos
1285 REAL(kind=dp), INTENT(IN) :: eps_lanczos
1286 REAL(kind=dp) :: homo, lumo
1287 INTEGER :: unit_nr
1288
1289 LOGICAL :: converged
1290 REAL(kind=dp) :: max_eig, min_eig, shift1, shift2
1291 TYPE(dbcsr_type) :: tmp1, tmp2, tmp3
1292
1293! temporary matrices used for HOMO/LUMO calculation
1294
1295 CALL dbcsr_create(tmp1, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1296
1297 CALL dbcsr_create(tmp2, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1298
1299 CALL dbcsr_create(tmp3, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1300
1301 shift1 = -eps_min
1302 shift2 = eps_max
1303
1304 ! find largest ev of P*(K+shift*1), where shift is the neg. val. of the smallest ev of K
1305 CALL dbcsr_copy(tmp2, matrix_k)
1306 CALL dbcsr_add_on_diag(tmp2, shift1)
1307 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, tmp2, &
1308 0.0_dp, tmp1, filter_eps=threshold)
1309 CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, &
1310 threshold=eps_lanczos, max_iter=max_iter_lanczos)
1311 homo = max_eig - shift1
1312 IF (unit_nr > 0) THEN
1313 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged
1314 END IF
1315
1316 ! -(1-P)*(K-shift*1) = (1-P)*(shift*1 - K), where shift is the largest ev of K
1317 CALL dbcsr_copy(tmp3, matrix_p)
1318 CALL dbcsr_scale(tmp3, -1.0_dp)
1319 CALL dbcsr_add_on_diag(tmp3, 1.0_dp) !tmp3 = 1-P
1320 CALL dbcsr_copy(tmp2, matrix_k)
1321 CALL dbcsr_add_on_diag(tmp2, -shift2)
1322 CALL dbcsr_multiply("N", "N", -1.0_dp, tmp3, tmp2, &
1323 0.0_dp, tmp1, filter_eps=threshold)
1324 CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, &
1325 threshold=eps_lanczos, max_iter=max_iter_lanczos)
1326 lumo = -max_eig + shift2
1327
1328 IF (unit_nr > 0) THEN
1329 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged
1330 WRITE (unit_nr, '(T6,A,1X,3F12.5)') 'HOMO/LUMO/gap', homo, lumo, lumo - homo
1331 END IF
1332 CALL dbcsr_release(tmp1)
1333 CALL dbcsr_release(tmp2)
1334 CALL dbcsr_release(tmp3)
1335
1336 END SUBROUTINE compute_homo_lumo
1337
1338! **************************************************************************************************
1339!> \brief ...
1340!> \param x ...
1341!> \param gamma_values ...
1342!> \param i ...
1343!> \return ...
1344! **************************************************************************************************
1345 FUNCTION evaluate_trs4_polynomial(x, gamma_values, i) RESULT(xr)
1346 REAL(kind=dp) :: x
1347 REAL(kind=dp), DIMENSION(:) :: gamma_values
1348 INTEGER :: i
1349 REAL(kind=dp) :: xr
1350
1351 REAL(kind=dp), PARAMETER :: gam_max = 6.0_dp, gam_min = 0.0_dp
1352
1353 INTEGER :: k
1354
1355 xr = x
1356 DO k = 1, i
1357 IF (gamma_values(k) > gam_max) THEN
1358 xr = 2*xr - xr**2
1359 ELSE IF (gamma_values(k) < gam_min) THEN
1360 xr = xr**2
1361 ELSE
1362 xr = (xr*xr)*(4*xr - 3*xr*xr) + gamma_values(k)*xr*xr*((1 - xr)**2)
1363 END IF
1364 END DO
1365 END FUNCTION evaluate_trs4_polynomial
1366
1367! **************************************************************************************************
1368!> \brief ...
1369!> \param homo ...
1370!> \param lumo ...
1371!> \param threshold ...
1372!> \return ...
1373! **************************************************************************************************
1374 FUNCTION estimate_steps(homo, lumo, threshold) RESULT(steps)
1375 REAL(kind=dp) :: homo, lumo, threshold
1376 INTEGER :: steps
1377
1378 INTEGER :: i
1379 REAL(kind=dp) :: h, l, m
1380
1381 l = lumo
1382 h = homo
1383
1384 DO i = 1, 200
1385 IF (abs(l) < threshold .AND. abs(1 - h) < threshold) EXIT
1386 m = 0.5_dp*(h + l)
1387 IF (m > 0.5_dp) THEN
1388 h = h**2
1389 l = l**2
1390 ELSE
1391 h = 2*h - h**2
1392 l = 2*l - l**2
1393 END IF
1394 END DO
1395 steps = i - 1
1396 END FUNCTION estimate_steps
1397
1398END MODULE dm_ls_scf_methods
arnoldi_api
arnoldi iteration using dbcsr
Definition arnoldi_api.F:15
arnoldi_api::arnoldi_extremal
subroutine, public arnoldi_extremal(matrix_a, max_ev, min_ev, converged, threshold, max_iter)
simple wrapper to estimate extremal eigenvalues with arnoldi, using the old lanczos interface this hi...
Definition arnoldi_api.F:254
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
dm_ls_scf_methods
lower level routines for linear scaling SCF
Definition dm_ls_scf_methods.F:14
dm_ls_scf_methods::ls_scf_init_matrix_s
subroutine, public ls_scf_init_matrix_s(matrix_s, ls_scf_env)
initialize S matrix related properties (sqrt, inverse...) Might be factored-out since this seems comm...
Definition dm_ls_scf_methods.F:72
dm_ls_scf_methods::density_matrix_sign_fixed_mu
subroutine, public density_matrix_sign_fixed_mu(matrix_p, trace, mu, sign_method, sign_order, matrix_ks, matrix_s, matrix_s_inv, threshold, sign_symmetric, submatrix_sign_method, matrix_s_sqrt_inv)
for a fixed mu, compute the corresponding density matrix and its trace
Definition dm_ls_scf_methods.F:484
dm_ls_scf_methods::apply_matrix_preconditioner
subroutine, public apply_matrix_preconditioner(matrix, direction, matrix_bs_sqrt, matrix_bs_sqrt_inv)
apply a preconditioner either forward (precondition) inv(sqrt(bs)) * A * inv(sqrt(bs)) backward (rest...
Definition dm_ls_scf_methods.F:316
dm_ls_scf_methods::density_matrix_sign
subroutine, public density_matrix_sign(matrix_p, mu, fixed_mu, sign_method, sign_order, matrix_ks, matrix_s, matrix_s_inv, nelectron, threshold, sign_symmetric, submatrix_sign_method, matrix_s_sqrt_inv)
compute the density matrix with a trace that is close to nelectron. take a mu as input,...
Definition dm_ls_scf_methods.F:374
dm_ls_scf_methods::compute_matrix_preconditioner
subroutine, public compute_matrix_preconditioner(matrix_s, preconditioner_type, ls_mstruct, matrix_bs_sqrt, matrix_bs_sqrt_inv, threshold, order, eps_lanczos, max_iter_lanczos)
compute for a block positive definite matrix s (bs) the sqrt(bs) and inv(sqrt(bs))
Definition dm_ls_scf_methods.F:215
dm_ls_scf_methods::compute_homo_lumo
subroutine, public compute_homo_lumo(matrix_k, matrix_p, eps_min, eps_max, threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr)
compute the homo and lumo given a KS matrix and a density matrix in the orthonormalized basis and the...
Definition dm_ls_scf_methods.F:1282
dm_ls_scf_methods::density_matrix_trs4
subroutine, public density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, e_mu, dynamic_threshold, matrix_ks_deviation, max_iter_lanczos, eps_lanczos, converged)
compute the density matrix using a trace-resetting algorithm
Definition dm_ls_scf_methods.F:723
dm_ls_scf_methods::density_matrix_tc2
subroutine, public density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, nelectron, threshold, e_homo, e_lumo, non_monotonic, eps_lanczos, max_iter_lanczos)
compute the density matrix using a non monotonic trace conserving algorithm based on SIAM DOI....
Definition dm_ls_scf_methods.F:1018
dm_ls_scf_qs
Routines for a linear scaling quickstep SCF run based on the density matrix, with a focus on the inte...
Definition dm_ls_scf_qs.F:15
dm_ls_scf_qs::matrix_qs_to_ls
subroutine, public matrix_qs_to_ls(matrix_ls, matrix_qs, ls_mstruct, covariant)
first link to QS, copy a QS matrix to LS matrix used to isolate QS style matrices from LS style will ...
Definition dm_ls_scf_qs.F:213
dm_ls_scf_types
Types needed for a linear scaling quickstep SCF run based on the density matrix.
Definition dm_ls_scf_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::ls_scf_submatrix_sign_direct_muadj
integer, parameter, public ls_scf_submatrix_sign_direct_muadj
Definition input_constants.F:965
input_constants::ls_s_preconditioner_molecular
integer, parameter, public ls_s_preconditioner_molecular
Definition input_constants.F:946
input_constants::ls_s_preconditioner_atomic
integer, parameter, public ls_s_preconditioner_atomic
Definition input_constants.F:946
input_constants::ls_scf_sign_submatrix
integer, parameter, public ls_scf_sign_submatrix
Definition input_constants.F:963
input_constants::ls_s_sqrt_proot
integer, parameter, public ls_s_sqrt_proot
Definition input_constants.F:957
input_constants::ls_s_sqrt_ns
integer, parameter, public ls_s_sqrt_ns
Definition input_constants.F:957
input_constants::ls_s_preconditioner_none
integer, parameter, public ls_s_preconditioner_none
Definition input_constants.F:946
input_constants::ls_scf_sign_proot
integer, parameter, public ls_scf_sign_proot
Definition input_constants.F:963
input_constants::ls_scf_submatrix_sign_ns
integer, parameter, public ls_scf_submatrix_sign_ns
Definition input_constants.F:965
input_constants::ls_scf_sign_ns
integer, parameter, public ls_scf_sign_ns
Definition input_constants.F:963
input_constants::ls_scf_submatrix_sign_direct_muadj_lowmem
integer, parameter, public ls_scf_submatrix_sign_direct_muadj_lowmem
Definition input_constants.F:965
input_constants::ls_cluster_atomic
integer, parameter, public ls_cluster_atomic
Definition input_constants.F:950
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::matrix_sqrt_newton_schulz
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Definition iterate_matrix.F:1624
iterate_matrix::matrix_sign_submatrix
subroutine, public matrix_sign_submatrix(matrix_sign, matrix, threshold, sign_order, submatrix_sign_method)
Submatrix method.
Definition iterate_matrix.F:1362
iterate_matrix::matrix_sign_submatrix_mu_adjust
subroutine, public matrix_sign_submatrix_mu_adjust(matrix_sign, matrix, mu, nelectron, threshold, variant)
Submatrix method with internal adjustment of chemical potential.
Definition iterate_matrix.F:1439
iterate_matrix::invert_hotelling
subroutine, public invert_hotelling(matrix_inverse, matrix, threshold, use_inv_as_guess, norm_convergence, filter_eps, accelerator_order, max_iter_lanczos, eps_lanczos, silent)
invert a symmetric positive definite matrix by Hotelling's method explicit symmetrization makes this ...
Definition iterate_matrix.F:472
iterate_matrix::matrix_sign_proot
subroutine, public matrix_sign_proot(matrix_sign, matrix, threshold, sign_order)
compute the sign a matrix using the general algorithm for the p-th root of Richters et al....
Definition iterate_matrix.F:1009
iterate_matrix::matrix_sign_newton_schulz
subroutine, public matrix_sign_newton_schulz(matrix_sign, matrix, threshold, sign_order)
compute the sign a matrix using Newton-Schulz iterations
Definition iterate_matrix.F:681
iterate_matrix::matrix_sqrt_proot
subroutine, public matrix_sqrt_proot(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
compute the sqrt of a matrix via the general algorithm for the p-th root of Richters et al....
Definition iterate_matrix.F:1871
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:106
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:123
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::abnormal_value
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Definition mathlib.F:150
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
dm_ls_scf_types::ls_mstruct_type
Definition dm_ls_scf_types.F:42
dm_ls_scf_types::ls_scf_env_type
Definition dm_ls_scf_types.F:88