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