(git:1f9fd2c)
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("")
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)) &
520 cpabort("Argument matrix_s_sqrt_inv required if sign_symmetric is set")
521
522 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
523 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv2, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
524 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
525 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, filter_eps=threshold)
526 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, matrix_s_sqrt_inv, &
527 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
528 CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)
529
530 SELECT CASE (sign_method)
531 CASE (ls_scf_sign_ns)
532 CALL matrix_sign_newton_schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, iounit=-1)
533 CASE (ls_scf_sign_proot)
534 CALL matrix_sign_proot(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
535 CASE (ls_scf_sign_submatrix)
536 CALL matrix_sign_submatrix(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, submatrix_sign_method)
537 CASE DEFAULT
538 cpabort("Unkown sign method.")
539 END SELECT
540 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
541 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2)
542
543 ELSE ! .NOT. sign_symmetric
544 ! get inv(S)*H-I*mu
545 CALL dbcsr_create(matrix_sinv_ks, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
546 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_inv, matrix_ks, &
547 0.0_dp, matrix_sinv_ks, filter_eps=threshold)
548 CALL dbcsr_add_on_diag(matrix_sinv_ks, -mu)
549
550 ! compute sign(inv(S)*H-I*mu)
551 SELECT CASE (sign_method)
552 CASE (ls_scf_sign_ns)
553 CALL matrix_sign_newton_schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order, iounit=-1)
554 CASE (ls_scf_sign_proot)
555 CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order)
556 CASE (ls_scf_sign_submatrix)
557 CALL matrix_sign_submatrix(matrix_sign, matrix_sinv_ks, threshold, sign_order, submatrix_sign_method)
558 CASE DEFAULT
559 cpabort("Unkown sign method.")
560 END SELECT
561 CALL dbcsr_release(matrix_sinv_ks)
562 END IF
563
564 ! now construct the density matrix PS=0.5*(I-sign(inv(S)H-I*mu))
565 CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
566 CALL dbcsr_copy(matrix_p_ud, matrix_sign)
567 CALL dbcsr_scale(matrix_p_ud, -0.5_dp)
568 CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp)
569 CALL dbcsr_release(matrix_sign)
570
571 ! we now have PS, lets get its trace
572 CALL dbcsr_trace(matrix_p_ud, trace)
573
574 ! we can also check it is idempotent PS*PS=PS
575 CALL dbcsr_create(matrix_tmp, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
576 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_p_ud, &
577 0.0_dp, matrix_tmp, filter_eps=threshold)
578 CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp)
579 frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
580 IF (unit_nr > 0 .AND. frob_matrix > 0.001_dp) &
581 WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix
582
583 IF (sign_symmetric) THEN
584 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, &
585 0.0_dp, matrix_tmp, filter_eps=threshold)
586 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_s_sqrt_inv, &
587 0.0_dp, matrix_p, filter_eps=threshold)
588 ELSE
589
590 ! get P=PS*inv(S)
591 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_s_inv, &
592 0.0_dp, matrix_p, filter_eps=threshold)
593 END IF
594 CALL dbcsr_release(matrix_p_ud)
595 CALL dbcsr_release(matrix_tmp)
596
597 CALL timestop(handle)
598
599 END SUBROUTINE density_matrix_sign_fixed_mu
600
601! **************************************************************************************************
602!> \brief compute the corresponding density matrix and its trace, using methods with internal mu adjustment
603!> \param matrix_p ...
604!> \param trace ...
605!> \param mu ...
606!> \param sign_method ...
607!> \param matrix_ks ...
608!> \param matrix_s ...
609!> \param threshold ...
610!> \param submatrix_sign_method ...
611!> \param nelectron ...
612!> \param matrix_s_sqrt_inv ...
613!> \par History
614!> 2020.07 created, based on density_matrix_sign_fixed_mu [Michael Lass]
615!> \author Michael Lass
616! **************************************************************************************************
617 SUBROUTINE density_matrix_sign_internal_mu(matrix_p, trace, mu, sign_method, matrix_ks, &
618 matrix_s, threshold, submatrix_sign_method, &
619 nelectron, matrix_s_sqrt_inv)
620
621 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
622 REAL(kind=dp), INTENT(OUT) :: trace
623 REAL(kind=dp), INTENT(INOUT) :: mu
624 INTEGER :: sign_method
625 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s
626 REAL(kind=dp), INTENT(IN) :: threshold
627 INTEGER :: submatrix_sign_method
628 INTEGER, INTENT(IN) :: nelectron
629 TYPE(dbcsr_type), INTENT(IN) :: matrix_s_sqrt_inv
630
631 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_sign_internal_mu'
632
633 INTEGER :: handle, unit_nr
634 REAL(kind=dp) :: frob_matrix
635 TYPE(cp_logger_type), POINTER :: logger
636 TYPE(dbcsr_type) :: matrix_p_ud, matrix_sign, &
637 matrix_ssqrtinv_ks_ssqrtinv, &
638 matrix_ssqrtinv_ks_ssqrtinv2, &
639 matrix_tmp
640
641 CALL timeset(routinen, handle)
642
643 logger => cp_get_default_logger()
644 IF (logger%para_env%is_source()) THEN
645 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
646 ELSE
647 unit_nr = -1
648 END IF
649
650 CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
651
652 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
653 CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv2, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
654 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
655 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, filter_eps=threshold)
656 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, matrix_s_sqrt_inv, &
657 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
658 CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)
659
660 SELECT CASE (sign_method)
661 CASE (ls_scf_sign_submatrix)
662 SELECT CASE (submatrix_sign_method)
663 CASE (ls_scf_submatrix_sign_direct_muadj, ls_scf_submatrix_sign_direct_muadj_lowmem)
664 CALL matrix_sign_submatrix_mu_adjust(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, mu, nelectron, threshold, &
665 submatrix_sign_method)
666 CASE DEFAULT
667 cpabort("density_matrix_sign_internal_mu called with invalid submatrix sign method")
668 END SELECT
669 CASE DEFAULT
670 cpabort("density_matrix_sign_internal_mu called with invalid sign method.")
671 END SELECT
672 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
673 CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2)
674
675 ! now construct the density matrix PS=0.5*(I-sign(inv(S)H-I*mu))
676 CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
677 CALL dbcsr_copy(matrix_p_ud, matrix_sign)
678 CALL dbcsr_scale(matrix_p_ud, -0.5_dp)
679 CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp)
680 CALL dbcsr_release(matrix_sign)
681
682 ! we now have PS, lets get its trace
683 CALL dbcsr_trace(matrix_p_ud, trace)
684
685 ! we can also check it is idempotent PS*PS=PS
686 CALL dbcsr_create(matrix_tmp, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
687 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_p_ud, &
688 0.0_dp, matrix_tmp, filter_eps=threshold)
689 CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp)
690 frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
691 IF (unit_nr > 0 .AND. frob_matrix > 0.001_dp) &
692 WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix
693
694 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, &
695 0.0_dp, matrix_tmp, filter_eps=threshold)
696 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_s_sqrt_inv, &
697 0.0_dp, matrix_p, filter_eps=threshold)
698 CALL dbcsr_release(matrix_p_ud)
699 CALL dbcsr_release(matrix_tmp)
700
701 CALL timestop(handle)
702
703 END SUBROUTINE density_matrix_sign_internal_mu
704
705! **************************************************************************************************
706!> \brief compute the density matrix using a trace-resetting algorithm
707!> \param matrix_p ...
708!> \param matrix_ks ...
709!> \param matrix_s_sqrt_inv ...
710!> \param nelectron ...
711!> \param threshold ...
712!> \param e_homo ...
713!> \param e_lumo ...
714!> \param e_mu ...
715!> \param dynamic_threshold ...
716!> \param matrix_ks_deviation ...
717!> \param max_iter_lanczos ...
718!> \param eps_lanczos ...
719!> \param converged ...
720!> \param iounit ...
721!> \par History
722!> 2012.06 created [Florian Thoele]
723!> \author Florian Thoele
724! **************************************************************************************************
725 SUBROUTINE density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, &
726 nelectron, threshold, e_homo, e_lumo, e_mu, &
727 dynamic_threshold, matrix_ks_deviation, &
728 max_iter_lanczos, eps_lanczos, converged, iounit)
729
730 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
731 TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv
732 INTEGER, INTENT(IN) :: nelectron
733 REAL(kind=dp), INTENT(IN) :: threshold
734 REAL(kind=dp), INTENT(INOUT) :: e_homo, e_lumo, e_mu
735 LOGICAL, INTENT(IN), OPTIONAL :: dynamic_threshold
736 TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: matrix_ks_deviation
737 INTEGER, INTENT(IN) :: max_iter_lanczos
738 REAL(kind=dp), INTENT(IN) :: eps_lanczos
739 LOGICAL, INTENT(OUT), OPTIONAL :: converged
740 INTEGER, INTENT(IN), OPTIONAL :: iounit
741
742 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_trs4'
743 INTEGER, PARAMETER :: max_iter = 100
744 REAL(kind=dp), PARAMETER :: gamma_max = 6.0_dp, gamma_min = 0.0_dp
745
746 INTEGER :: branch, estimated_steps, handle, i, j, &
747 unit_nr
748 INTEGER(kind=int_8) :: flop1, flop2
749 LOGICAL :: arnoldi_converged, do_dyn_threshold
750 REAL(kind=dp) :: current_threshold, delta_n, eps_max, eps_min, est_threshold, frob_id, &
751 frob_x, gam, homo, lumo, max_eig, max_threshold, maxdev, maxev, min_eig, minev, mmin, mu, &
752 mu_a, mu_b, mu_c, mu_fa, mu_fc, occ_matrix, scaled_homo_bound, scaled_lumo_bound, t1, t2, &
753 trace_fx, trace_gx, xi
754 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: gamma_values
755 TYPE(cp_logger_type), POINTER :: logger
756 TYPE(dbcsr_type) :: matrix_k0, matrix_x, matrix_x_nosym, &
757 matrix_xidsq, matrix_xsq, tmp_gx
758
759 IF (nelectron == 0) THEN
760 CALL dbcsr_set(matrix_p, 0.0_dp)
761 RETURN
762 END IF
763
764 CALL timeset(routinen, handle)
765
766 IF (PRESENT(iounit)) THEN
767 unit_nr = iounit
768 ELSE
769 logger => cp_get_default_logger()
770 IF (logger%para_env%is_source()) THEN
771 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
772 ELSE
773 unit_nr = -1
774 END IF
775 END IF
776
777 do_dyn_threshold = .false.
778 IF (PRESENT(dynamic_threshold)) do_dyn_threshold = dynamic_threshold
779
780 IF (PRESENT(converged)) converged = .false.
781
782 ! init X = (eps_n*I - H)/(eps_n - eps_0) ... H* = S^-1/2*H*S^-1/2
783 CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type="S")
784
785 ! at some points the non-symmetric version of x is required
786 CALL dbcsr_create(matrix_x_nosym, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
787
788 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
789 0.0_dp, matrix_x_nosym, filter_eps=threshold)
790 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x_nosym, matrix_s_sqrt_inv, &
791 0.0_dp, matrix_x, filter_eps=threshold)
792 CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym)
793
794 CALL dbcsr_create(matrix_k0, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
795 CALL dbcsr_copy(matrix_k0, matrix_x_nosym)
796
797 ! compute the deviation in the mixed matrix, as seen in the ortho basis
798 IF (do_dyn_threshold) THEN
799 cpassert(PRESENT(matrix_ks_deviation))
800 CALL dbcsr_add(matrix_ks_deviation, matrix_x_nosym, -1.0_dp, 1.0_dp)
801 CALL arnoldi_extremal(matrix_ks_deviation, maxev, minev, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
802 converged=arnoldi_converged)
803 maxdev = max(abs(maxev), abs(minev))
804 IF (unit_nr > 0) THEN
805 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", arnoldi_converged
806 WRITE (unit_nr, '(T6,A,1X,F12.5)') "change in mixed matrix: ", maxdev
807 WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo + maxdev
808 WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo - maxdev
809 WRITE (unit_nr, '(T6,A,1X,L12)') "Predicts a gap ? ", ((e_lumo - maxdev) - (e_homo + maxdev)) > 0
810 END IF
811 ! save the old mixed matrix
812 CALL dbcsr_copy(matrix_ks_deviation, matrix_x_nosym)
813
814 END IF
815
816 ! get largest/smallest eigenvalues for scaling
817 CALL arnoldi_extremal(matrix_x_nosym, max_eig, min_eig, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
818 converged=arnoldi_converged)
819 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", &
820 min_eig, max_eig, " converged: ", arnoldi_converged
821 eps_max = max_eig
822 eps_min = min_eig
823
824 ! scale KS matrix
825 IF (eps_max == eps_min) THEN
826 CALL dbcsr_scale(matrix_x, 1.0_dp/eps_max)
827 ELSE
828 CALL dbcsr_add_on_diag(matrix_x, -eps_max)
829 CALL dbcsr_scale(matrix_x, -1.0_dp/(eps_max - eps_min))
830 END IF
831
832 current_threshold = threshold
833 IF (do_dyn_threshold) THEN
834 ! scale bounds for HOMO/LUMO
835 scaled_homo_bound = (eps_max - (e_homo + maxdev))/(eps_max - eps_min)
836 scaled_lumo_bound = (eps_max - (e_lumo - maxdev))/(eps_max - eps_min)
837 END IF
838
839 CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type="S")
840
841 CALL dbcsr_create(matrix_xidsq, template=matrix_ks, matrix_type="S")
842
843 CALL dbcsr_create(tmp_gx, template=matrix_ks, matrix_type="S")
844
845 ALLOCATE (gamma_values(max_iter))
846
847 DO i = 1, max_iter
848 t1 = m_walltime()
849 flop1 = 0; flop2 = 0
850
851 ! get X*X
852 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, &
853 0.0_dp, matrix_xsq, &
854 filter_eps=current_threshold, flop=flop1)
855
856 ! intermediate use matrix_xidsq to compute = X*X-X
857 CALL dbcsr_copy(matrix_xidsq, matrix_x)
858 CALL dbcsr_add(matrix_xidsq, matrix_xsq, -1.0_dp, 1.0_dp)
859 frob_id = dbcsr_frobenius_norm(matrix_xidsq)
860 frob_x = dbcsr_frobenius_norm(matrix_x)
861
862 ! xidsq = (1-X)*(1-X)
863 ! use (1-x)*(1-x) = 1 + x*x - 2*x
864 CALL dbcsr_copy(matrix_xidsq, matrix_x)
865 CALL dbcsr_add(matrix_xidsq, matrix_xsq, -2.0_dp, 1.0_dp)
866 CALL dbcsr_add_on_diag(matrix_xidsq, 1.0_dp)
867
868 ! tmp_gx = 4X-3X*X
869 CALL dbcsr_copy(tmp_gx, matrix_x)
870 CALL dbcsr_add(tmp_gx, matrix_xsq, 4.0_dp, -3.0_dp)
871
872 ! get gamma
873 ! Tr(F) = Tr(XX*tmp_gx) Tr(G) is equivalent
874 CALL dbcsr_dot(matrix_xsq, matrix_xidsq, trace_gx)
875 CALL dbcsr_dot(matrix_xsq, tmp_gx, trace_fx)
876
877 ! if converged, and gam becomes noisy, fix it to 3, which results in a final McWeeny step.
878 ! do this only if the electron count is reasonable.
879 ! maybe tune if the current criterion is not good enough
880 delta_n = nelectron - trace_fx
881 ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ...
882 IF (((frob_id*frob_id) < (threshold*frob_x*frob_x)) .AND. (abs(delta_n) < 0.5_dp)) THEN
883 gam = 3.0_dp
884 ELSE IF (abs(delta_n) < 1e-14_dp) THEN
885 gam = 0.0_dp ! rare case of perfect electron count
886 ELSE
887 ! make sure, we don't divide by zero, as soon as gam is outside the interval gam_min,gam_max, it doesn't matter
888 gam = delta_n/max(trace_gx, abs(delta_n)/100)
889 END IF
890 gamma_values(i) = gam
891
892 IF (unit_nr > 0 .AND. .false.) THEN
893 WRITE (unit_nr, *) "trace_fx", trace_fx, "trace_gx", trace_gx, "gam", gam, &
894 "frob_id", frob_id, "conv", abs(frob_id/frob_x)
895 END IF
896
897 IF (do_dyn_threshold) THEN
898 ! quantities used for dynamic thresholding, when the estimated gap is larger than zero
899 xi = (scaled_homo_bound - scaled_lumo_bound)
900 IF (xi > 0.0_dp) THEN
901 mmin = 0.5*(scaled_homo_bound + scaled_lumo_bound)
902 max_threshold = abs(1 - 2*mmin)*xi
903
904 scaled_homo_bound = evaluate_trs4_polynomial(scaled_homo_bound, gamma_values(i:), 1)
905 scaled_lumo_bound = evaluate_trs4_polynomial(scaled_lumo_bound, gamma_values(i:), 1)
906 estimated_steps = estimate_steps(scaled_homo_bound, scaled_lumo_bound, threshold)
907
908 est_threshold = (threshold/(estimated_steps + i + 1))*xi/(1 + threshold/(estimated_steps + i + 1))
909 est_threshold = min(max_threshold, est_threshold)
910 IF (i > 1) est_threshold = max(est_threshold, 0.1_dp*current_threshold)
911 current_threshold = est_threshold
912 ELSE
913 current_threshold = threshold
914 END IF
915 END IF
916
917 IF (gam > gamma_max) THEN
918 ! Xn+1 = 2X-X*X
919 CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp)
920 CALL dbcsr_filter(matrix_x, current_threshold)
921 branch = 1
922 ELSE IF (gam < gamma_min) THEN
923 ! Xn+1 = X*X
924 CALL dbcsr_copy(matrix_x, matrix_xsq)
925 branch = 2
926 ELSE
927 ! Xn+1 = F(X) + gam*G(X)
928 CALL dbcsr_add(tmp_gx, matrix_xidsq, 1.0_dp, gam)
929 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_xsq, tmp_gx, &
930 0.0_dp, matrix_x, &
931 flop=flop2, filter_eps=current_threshold)
932 branch = 3
933 END IF
934
935 occ_matrix = dbcsr_get_occupation(matrix_x)
936 t2 = m_walltime()
937 IF (unit_nr > 0) THEN
938 WRITE (unit_nr, &
939 '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TRS4 it ", &
940 i, occ_matrix, abs(trace_gx), t2 - t1, &
941 (flop1 + flop2)/(1.0e6_dp*max(t2 - t1, 0.001_dp)), current_threshold
942 CALL m_flush(unit_nr)
943 END IF
944
945 IF (abnormal_value(trace_gx)) &
946 cpabort("trace_gx is an abnormal value (NaN/Inf).")
947
948 ! a branch of 1 or 2 appears to lead to a less accurate electron number count and premature exit
949 ! if it turns out this does not exit because we get stuck in branch 1/2 for a reason we need to refine further
950 ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ...
951 IF ((frob_id*frob_id) < (threshold*frob_x*frob_x) .AND. branch == 3 .AND. (abs(delta_n) < 0.5_dp)) THEN
952 IF (PRESENT(converged)) converged = .true.
953 EXIT
954 END IF
955
956 END DO
957
958 occ_matrix = dbcsr_get_occupation(matrix_x)
959 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,I3,1X,F10.8,E12.3)') 'Final TRS4 iteration ', i, occ_matrix, abs(trace_gx)
960
961 ! free some memory
962 CALL dbcsr_release(tmp_gx)
963 CALL dbcsr_release(matrix_xsq)
964 CALL dbcsr_release(matrix_xidsq)
965
966 ! output to matrix_p, P = inv(S)^0.5 X inv(S)^0.5
967 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, &
968 0.0_dp, matrix_x_nosym, filter_eps=threshold)
969 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_x_nosym, &
970 0.0_dp, matrix_p, filter_eps=threshold)
971
972 ! calculate the chemical potential by doing a bisection of fk(x0)-0.5, where fk is evaluated using the stored values for gamma
973 ! E. Rubensson et al., Chem Phys Lett 432, 2006, 591-594
974 mu_a = 0.0_dp; mu_b = 1.0_dp;
975 mu_fa = evaluate_trs4_polynomial(mu_a, gamma_values, i - 1) - 0.5_dp
976 DO j = 1, 40
977 mu_c = 0.5*(mu_a + mu_b)
978 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
979 IF (abs(mu_fc) < 1.0e-6_dp .OR. (mu_b - mu_a)/2 < 1.0e-6_dp) EXIT !TODO: define threshold values
980
981 IF (mu_fc*mu_fa > 0) THEN
982 mu_a = mu_c
983 mu_fa = mu_fc
984 ELSE
985 mu_b = mu_c
986 END IF
987 END DO
988 mu = (eps_min - eps_max)*mu_c + eps_max
989 DEALLOCATE (gamma_values)
990 IF (unit_nr > 0) THEN
991 WRITE (unit_nr, '(T6,A,1X,F12.5)') 'Chemical potential (mu): ', mu
992 END IF
993 e_mu = mu
994
995 IF (do_dyn_threshold) THEN
996 CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym)
997 CALL compute_homo_lumo(matrix_k0, matrix_x_nosym, eps_min, eps_max, &
998 threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr)
999 e_homo = homo
1000 e_lumo = lumo
1001 END IF
1002
1003 CALL dbcsr_release(matrix_x)
1004 CALL dbcsr_release(matrix_x_nosym)
1005 CALL dbcsr_release(matrix_k0)
1006 CALL timestop(handle)
1007
1008 END SUBROUTINE density_matrix_trs4
1009
1010! **************************************************************************************************
1011!> \brief compute the density matrix using a non monotonic trace conserving
1012!> algorithm based on SIAM DOI. 10.1137/130911585.
1013!> 2014.04 created [Jonathan Mullin]
1014!> \param matrix_p ...
1015!> \param matrix_ks ...
1016!> \param matrix_s_sqrt_inv ...
1017!> \param nelectron ...
1018!> \param threshold ...
1019!> \param e_homo ...
1020!> \param e_lumo ...
1021!> \param non_monotonic ...
1022!> \param eps_lanczos ...
1023!> \param max_iter_lanczos ...
1024!> \param iounit ...
1025!> \author Jonathan Mullin
1026! **************************************************************************************************
1027 SUBROUTINE density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, &
1028 nelectron, threshold, e_homo, e_lumo, &
1029 non_monotonic, eps_lanczos, max_iter_lanczos, iounit)
1030
1031 TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
1032 TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv
1033 INTEGER, INTENT(IN) :: nelectron
1034 REAL(kind=dp), INTENT(IN) :: threshold
1035 REAL(kind=dp), INTENT(INOUT) :: e_homo, e_lumo
1036 LOGICAL, INTENT(IN), OPTIONAL :: non_monotonic
1037 REAL(kind=dp), INTENT(IN) :: eps_lanczos
1038 INTEGER, INTENT(IN) :: max_iter_lanczos
1039 INTEGER, INTENT(IN), OPTIONAL :: iounit
1040
1041 CHARACTER(LEN=*), PARAMETER :: routinen = 'density_matrix_tc2'
1042 INTEGER, PARAMETER :: max_iter = 100
1043
1044 INTEGER :: handle, i, j, k, unit_nr
1045 INTEGER(kind=int_8) :: flop1, flop2
1046 LOGICAL :: converged, do_non_monotonic, &
1047 tc2_converged
1048 REAL(kind=dp) :: beta, betab, eps_max, eps_min, gama, &
1049 max_eig, min_eig, occ_matrix, t1, t2, &
1050 trace_fx, trace_gx
1051 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: alpha, lambda, nu, poly, wu, x, y
1052 TYPE(cp_logger_type), POINTER :: logger
1053 TYPE(dbcsr_type) :: matrix_tmp, matrix_x, matrix_xsq
1054
1055 CALL timeset(routinen, handle)
1056
1057 IF (PRESENT(iounit)) THEN
1058 unit_nr = iounit
1059 ELSE
1060 logger => cp_get_default_logger()
1061 IF (logger%para_env%is_source()) THEN
1062 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
1063 ELSE
1064 unit_nr = -1
1065 END IF
1066 END IF
1067
1068 do_non_monotonic = .false.
1069 IF (PRESENT(non_monotonic)) do_non_monotonic = non_monotonic
1070
1071 ! init X = (eps_n*I - H)/(eps_n - eps_0) ... H* = S^-1/2*H*S^-1/2
1072 CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1073 CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1074
1075 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
1076 0.0_dp, matrix_xsq, filter_eps=threshold)
1077 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_xsq, matrix_s_sqrt_inv, &
1078 0.0_dp, matrix_x, filter_eps=threshold)
1079
1080 IF (unit_nr > 0) THEN
1081 WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo
1082 WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo
1083 WRITE (unit_nr, '(T6,A,1X,L12)') "Predicts a gap ? ", ((e_lumo) - (e_homo)) > 0
1084 END IF
1085
1086 ! get largest/smallest eigenvalues for scaling
1087 CALL arnoldi_extremal(matrix_x, max_eig, min_eig, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
1088 converged=converged)
1089 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", &
1090 min_eig, max_eig, " converged: ", converged
1091
1092 eps_max = max_eig
1093 eps_min = min_eig
1094
1095 ! scale KS matrix
1096 CALL dbcsr_scale(matrix_x, -1.0_dp)
1097 CALL dbcsr_add_on_diag(matrix_x, eps_max)
1098 CALL dbcsr_scale(matrix_x, 1/(eps_max - eps_min))
1099
1100 CALL dbcsr_copy(matrix_xsq, matrix_x)
1101
1102 CALL dbcsr_create(matrix_tmp, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1103
1104 ALLOCATE (poly(max_iter))
1105 ALLOCATE (nu(max_iter))
1106 ALLOCATE (wu(max_iter))
1107 ALLOCATE (alpha(max_iter))
1108 ALLOCATE (x(4))
1109 ALLOCATE (y(4))
1110 ALLOCATE (lambda(4))
1111
1112! Controls over the non_monotonic bounds, First if low gap, bias slightly
1113 beta = (eps_max - abs(e_lumo))/(eps_max - eps_min)
1114 betab = (eps_max + abs(e_homo))/(eps_max - eps_min)
1115
1116 IF ((beta - betab) < 0.005_dp) THEN
1117 beta = beta - 0.002_dp
1118 betab = betab + 0.002_dp
1119 END IF
1120! Check if input specifies to use monotonic bounds.
1121 IF (.NOT. do_non_monotonic) THEN
1122 beta = 0.0_dp
1123 betab = 1.0_dp
1124 END IF
1125! initial SCF cycle has no reliable estimate of homo/lumo, force monotinic bounds.
1126 IF (e_homo == 0.0_dp) THEN
1127 beta = 0.0_dp
1128 betab = 1.0_dp
1129 END IF
1130
1131 ! init to take true branch first
1132 trace_fx = nelectron
1133 trace_gx = 0
1134
1135 tc2_converged = .false.
1136 DO i = 1, max_iter
1137 t1 = m_walltime()
1138 flop1 = 0; flop2 = 0
1139
1140 IF (abs(trace_fx - nelectron) <= abs(trace_gx - nelectron)) THEN
1141! Xn+1 = (aX+ (1-a)I)^2
1142 poly(i) = 1.0_dp
1143 alpha(i) = 2.0_dp/(2.0_dp - beta)
1144
1145 CALL dbcsr_scale(matrix_x, alpha(i))
1146 CALL dbcsr_add_on_diag(matrix_x, 1.0_dp - 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_copy(matrix_x, matrix_xsq)
1155
1156 beta = (1.0_dp - alpha(i)) + alpha(i)*beta
1157 beta = beta*beta
1158 betab = (1.0_dp - alpha(i)) + alpha(i)*betab
1159 betab = betab*betab
1160 ELSE
1161! Xn+1 = 2aX-a^2*X^2
1162 poly(i) = 0.0_dp
1163 alpha(i) = 2.0_dp/(1.0_dp + betab)
1164
1165 CALL dbcsr_scale(matrix_x, alpha(i))
1166 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, &
1167 0.0_dp, matrix_xsq, &
1168 filter_eps=threshold, flop=flop1)
1169
1170!save X for control variables
1171 CALL dbcsr_copy(matrix_tmp, matrix_x)
1172!
1173 CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp)
1174
1175 beta = alpha(i)*beta
1176 beta = 2.0_dp*beta - beta*beta
1177 betab = alpha(i)*betab
1178 betab = 2.0_dp*betab - betab*betab
1179
1180 END IF
1181 occ_matrix = dbcsr_get_occupation(matrix_x)
1182 t2 = m_walltime()
1183 IF (unit_nr > 0) THEN
1184 WRITE (unit_nr, &
1185 '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TC2 it ", &
1186 i, occ_matrix, t2 - t1, &
1187 (flop1 + flop2)/(1.0e6_dp*(t2 - t1)), threshold
1188 CALL m_flush(unit_nr)
1189 END IF
1190
1191! calculate control terms
1192 CALL dbcsr_trace(matrix_xsq, trace_fx)
1193
1194! intermediate use matrix_xsq compute X- X*X , temorarily use trace_gx
1195 CALL dbcsr_add(matrix_xsq, matrix_tmp, -1.0_dp, 1.0_dp)
1196 CALL dbcsr_trace(matrix_xsq, trace_gx)
1197 nu(i) = dbcsr_frobenius_norm(matrix_xsq)
1198 wu(i) = trace_gx
1199
1200! intermediate use matrix_xsq to compute = 2X - X*X
1201 CALL dbcsr_add(matrix_xsq, matrix_tmp, 1.0_dp, 1.0_dp)
1202 CALL dbcsr_trace(matrix_xsq, trace_gx)
1203! TC2 has quadratic convergence, using the frobeniums norm as an idempotency deviation test.
1204 IF (abs(nu(i)) < (threshold)) THEN
1205 tc2_converged = .true.
1206 EXIT
1207 END IF
1208 END DO
1209 IF (.NOT. tc2_converged) i = max_iter
1210
1211 occ_matrix = dbcsr_get_occupation(matrix_x)
1212 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))
1213
1214 ! output to matrix_p, P = inv(S)^0.5 X inv(S)^0.5
1215 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, &
1216 0.0_dp, matrix_tmp, filter_eps=threshold)
1217 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_tmp, &
1218 0.0_dp, matrix_p, filter_eps=threshold)
1219
1220 CALL dbcsr_release(matrix_xsq)
1221 CALL dbcsr_release(matrix_tmp)
1222
1223 ! ALGO 3 from. SIAM DOI. 10.1137/130911585
1224 x(1) = 1.0_dp
1225 x(2) = 1.0_dp
1226 x(3) = 0.0_dp
1227 x(4) = 0.0_dp
1228 gama = 6.0_dp - 4.0_dp*(sqrt(2.0_dp))
1229 gama = gama - gama*gama
1230 DO
1231 IF (i < 1) EXIT
1232 IF (nu(i) >= gama) EXIT
1233 ! safeguard against negative root, is skipping correct?
1234 IF (wu(i) < 1.0e-14_dp) THEN
1235 i = i - 1
1236 cycle
1237 END IF
1238 IF ((1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)) < 0.0_dp) THEN
1239 i = i - 1
1240 cycle
1241 END IF
1242 y(1) = 0.5_dp*(1.0_dp - sqrt(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)))
1243 y(2) = 0.5_dp*(1.0_dp - sqrt(1.0_dp - 4.0_dp*nu(i)))
1244 y(3) = 0.5_dp*(1.0_dp + sqrt(1.0_dp - 4.0_dp*nu(i)))
1245 y(4) = 0.5_dp*(1.0_dp + sqrt(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)))
1246 y(:) = min(1.0_dp, max(0.0_dp, y(:)))
1247 DO j = i, 1, -1
1248 IF (poly(j) == 1.0_dp) THEN
1249 DO k = 1, 4
1250 y(k) = sqrt(y(k))
1251 y(k) = (y(k) - 1.0_dp + alpha(j))/alpha(j)
1252 END DO ! end K
1253 ELSE
1254 DO k = 1, 4
1255 y(k) = 1.0_dp - sqrt(1.0_dp - y(k))
1256 y(k) = y(k)/alpha(j)
1257 END DO ! end K
1258 END IF ! end poly
1259 END DO ! end j
1260 x(1) = min(x(1), y(1))
1261 x(2) = min(x(2), y(2))
1262 x(3) = max(x(3), y(3))
1263 x(4) = max(x(4), y(4))
1264 i = i - 1
1265 END DO ! end i
1266! lambda 1,2,3,4 are:: out lumo, in lumo, in homo, out homo
1267 DO k = 1, 4
1268 lambda(k) = eps_max - (eps_max - eps_min)*x(k)
1269 END DO ! end k
1270! END ALGO 3 from. SIAM DOI. 10.1137/130911585
1271 e_homo = lambda(4)
1272 e_lumo = lambda(1)
1273 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "outer homo/lumo/gap", e_homo, e_lumo, (e_lumo - e_homo)
1274 IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "inner homo/lumo/gap", lambda(3), lambda(2), (lambda(2) - lambda(3))
1275
1276 DEALLOCATE (poly)
1277 DEALLOCATE (nu)
1278 DEALLOCATE (wu)
1279 DEALLOCATE (alpha)
1280 DEALLOCATE (x)
1281 DEALLOCATE (y)
1282 DEALLOCATE (lambda)
1283
1284 CALL dbcsr_release(matrix_x)
1285 CALL timestop(handle)
1286
1287 END SUBROUTINE density_matrix_tc2
1288
1289! **************************************************************************************************
1290!> \brief compute the homo and lumo given a KS matrix and a density matrix in the orthonormalized basis
1291!> and the eps_min and eps_max, min and max eigenvalue of the ks matrix
1292!> \param matrix_k ...
1293!> \param matrix_p ...
1294!> \param eps_min ...
1295!> \param eps_max ...
1296!> \param threshold ...
1297!> \param max_iter_lanczos ...
1298!> \param eps_lanczos ...
1299!> \param homo ...
1300!> \param lumo ...
1301!> \param unit_nr ...
1302!> \par History
1303!> 2012.06 created [Florian Thoele]
1304!> \author Florian Thoele
1305! **************************************************************************************************
1306 SUBROUTINE compute_homo_lumo(matrix_k, matrix_p, eps_min, eps_max, threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr)
1307 TYPE(dbcsr_type) :: matrix_k, matrix_p
1308 REAL(kind=dp) :: eps_min, eps_max, threshold
1309 INTEGER, INTENT(IN) :: max_iter_lanczos
1310 REAL(kind=dp), INTENT(IN) :: eps_lanczos
1311 REAL(kind=dp) :: homo, lumo
1312 INTEGER :: unit_nr
1313
1314 LOGICAL :: converged
1315 REAL(kind=dp) :: max_eig, min_eig, shift1, shift2
1316 TYPE(dbcsr_type) :: tmp1, tmp2, tmp3
1317
1318! temporary matrices used for HOMO/LUMO calculation
1319
1320 CALL dbcsr_create(tmp1, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1321
1322 CALL dbcsr_create(tmp2, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1323
1324 CALL dbcsr_create(tmp3, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1325
1326 shift1 = -eps_min
1327 shift2 = eps_max
1328
1329 ! find largest ev of P*(K+shift*1), where shift is the neg. val. of the smallest ev of K
1330 CALL dbcsr_copy(tmp2, matrix_k)
1331 CALL dbcsr_add_on_diag(tmp2, shift1)
1332 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, tmp2, &
1333 0.0_dp, tmp1, filter_eps=threshold)
1334 CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, &
1335 threshold=eps_lanczos, max_iter=max_iter_lanczos)
1336 homo = max_eig - shift1
1337 IF (unit_nr > 0) THEN
1338 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged
1339 END IF
1340
1341 ! -(1-P)*(K-shift*1) = (1-P)*(shift*1 - K), where shift is the largest ev of K
1342 CALL dbcsr_copy(tmp3, matrix_p)
1343 CALL dbcsr_scale(tmp3, -1.0_dp)
1344 CALL dbcsr_add_on_diag(tmp3, 1.0_dp) !tmp3 = 1-P
1345 CALL dbcsr_copy(tmp2, matrix_k)
1346 CALL dbcsr_add_on_diag(tmp2, -shift2)
1347 CALL dbcsr_multiply("N", "N", -1.0_dp, tmp3, tmp2, &
1348 0.0_dp, tmp1, filter_eps=threshold)
1349 CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, &
1350 threshold=eps_lanczos, max_iter=max_iter_lanczos)
1351 lumo = -max_eig + shift2
1352
1353 IF (unit_nr > 0) THEN
1354 WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged
1355 WRITE (unit_nr, '(T6,A,1X,3F12.5)') 'HOMO/LUMO/gap', homo, lumo, lumo - homo
1356 END IF
1357 CALL dbcsr_release(tmp1)
1358 CALL dbcsr_release(tmp2)
1359 CALL dbcsr_release(tmp3)
1360
1361 END SUBROUTINE compute_homo_lumo
1362
1363! **************************************************************************************************
1364!> \brief ...
1365!> \param x ...
1366!> \param gamma_values ...
1367!> \param i ...
1368!> \return ...
1369! **************************************************************************************************
1370 FUNCTION evaluate_trs4_polynomial(x, gamma_values, i) RESULT(xr)
1371 REAL(kind=dp) :: x
1372 REAL(kind=dp), DIMENSION(:) :: gamma_values
1373 INTEGER :: i
1374 REAL(kind=dp) :: xr
1375
1376 REAL(kind=dp), PARAMETER :: gam_max = 6.0_dp, gam_min = 0.0_dp
1377
1378 INTEGER :: k
1379
1380 xr = x
1381 DO k = 1, i
1382 IF (gamma_values(k) > gam_max) THEN
1383 xr = 2*xr - xr**2
1384 ELSE IF (gamma_values(k) < gam_min) THEN
1385 xr = xr**2
1386 ELSE
1387 xr = (xr*xr)*(4*xr - 3*xr*xr) + gamma_values(k)*xr*xr*((1 - xr)**2)
1388 END IF
1389 END DO
1390 END FUNCTION evaluate_trs4_polynomial
1391
1392! **************************************************************************************************
1393!> \brief ...
1394!> \param homo ...
1395!> \param lumo ...
1396!> \param threshold ...
1397!> \return ...
1398! **************************************************************************************************
1399 FUNCTION estimate_steps(homo, lumo, threshold) RESULT(steps)
1400 REAL(kind=dp) :: homo, lumo, threshold
1401 INTEGER :: steps
1402
1403 INTEGER :: i
1404 REAL(kind=dp) :: h, l, m
1405
1406 l = lumo
1407 h = homo
1408
1409 DO i = 1, 200
1410 IF (abs(l) < threshold .AND. abs(1 - h) < threshold) EXIT
1411 m = 0.5_dp*(h + l)
1412 IF (m > 0.5_dp) THEN
1413 h = h**2
1414 l = l**2
1415 ELSE
1416 h = 2*h - h**2
1417 l = 2*l - l**2
1418 END IF
1419 END DO
1420 steps = i - 1
1421 END FUNCTION estimate_steps
1422
1423END 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:1187
cp_dbcsr_api::dbcsr_iterator_blocks_left
logical function, public dbcsr_iterator_blocks_left(iterator)
...
Definition cp_dbcsr_api.F:943
cp_dbcsr_api::dbcsr_iterator_stop
subroutine, public dbcsr_iterator_stop(iterator)
...
Definition cp_dbcsr_api.F:1049
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:524
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:370
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:1095
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:970
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:657
cp_dbcsr_api::dbcsr_get_occupation
real(kind=dp) function, public dbcsr_get_occupation(matrix)
...
Definition cp_dbcsr_api.F:892
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:672
cp_dbcsr_api::dbcsr_iterator_start
subroutine, public dbcsr_iterator_start(iterator, matrix, shared, dynamic, dynamic_byrows)
...
Definition cp_dbcsr_api.F:1011
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1203
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1141
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1124
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:253
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: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::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:729
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:1030
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:1307
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:1013
input_constants::ls_s_preconditioner_molecular
integer, parameter, public ls_s_preconditioner_molecular
Definition input_constants.F:994
input_constants::ls_s_preconditioner_atomic
integer, parameter, public ls_s_preconditioner_atomic
Definition input_constants.F:994
input_constants::ls_scf_sign_submatrix
integer, parameter, public ls_scf_sign_submatrix
Definition input_constants.F:1011
input_constants::ls_s_sqrt_proot
integer, parameter, public ls_s_sqrt_proot
Definition input_constants.F:1005
input_constants::ls_s_sqrt_ns
integer, parameter, public ls_s_sqrt_ns
Definition input_constants.F:1005
input_constants::ls_s_preconditioner_none
integer, parameter, public ls_s_preconditioner_none
Definition input_constants.F:994
input_constants::ls_scf_sign_proot
integer, parameter, public ls_scf_sign_proot
Definition input_constants.F:1011
input_constants::ls_scf_submatrix_sign_ns
integer, parameter, public ls_scf_submatrix_sign_ns
Definition input_constants.F:1013
input_constants::ls_scf_sign_ns
integer, parameter, public ls_scf_sign_ns
Definition input_constants.F:1011
input_constants::ls_scf_submatrix_sign_direct_muadj_lowmem
integer, parameter, public ls_scf_submatrix_sign_direct_muadj_lowmem
Definition input_constants.F:1013
input_constants::ls_cluster_atomic
integer, parameter, public ls_cluster_atomic
Definition input_constants.F:998
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:1365
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:1442
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:1626
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:1878
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:158
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