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