Line data Source code
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 : ! **************************************************************************************************
14 : MODULE 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 :
62 : CONTAINS
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 12822 : 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 12822 : CALL timeset(routineN, handle)
85 :
86 : ! get a useful output_unit
87 12822 : logger => cp_get_default_logger()
88 12822 : IF (logger%para_env%is_source()) THEN
89 6411 : 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 12822 : IF (ls_scf_env%has_unit_metric) THEN
96 14 : CALL dbcsr_set(ls_scf_env%matrix_s, 0.0_dp)
97 14 : CALL dbcsr_add_on_diag(ls_scf_env%matrix_s, 1.0_dp)
98 : ELSE
99 12808 : CALL matrix_qs_to_ls(ls_scf_env%matrix_s, matrix_s, ls_scf_env%ls_mstruct, covariant=.TRUE.)
100 : END IF
101 :
102 12822 : CALL dbcsr_filter(ls_scf_env%matrix_s, ls_scf_env%eps_filter)
103 :
104 : ! needs a preconditioner for S
105 12822 : 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 456 : matrix_type=dbcsr_type_no_symmetry)
108 : CALL dbcsr_create(ls_scf_env%matrix_bs_sqrt_inv, template=ls_scf_env%matrix_s, &
109 456 : 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 456 : ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos)
115 : END IF
116 :
117 : ! precondition S
118 12822 : IF (ls_scf_env%has_s_preconditioner) THEN
119 : CALL apply_matrix_preconditioner(ls_scf_env%matrix_s, "forward", &
120 456 : 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 12822 : 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 12818 : matrix_type=dbcsr_type_no_symmetry)
128 : CALL dbcsr_create(ls_scf_env%matrix_s_sqrt_inv, template=ls_scf_env%matrix_s, &
129 12818 : matrix_type=dbcsr_type_no_symmetry)
130 :
131 12826 : 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 8 : 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 12810 : iounit=-1)
144 : CASE DEFAULT
145 12818 : CPABORT("Unknown sqrt method.")
146 : END SELECT
147 :
148 12818 : IF (ls_scf_env%check_s_inv) THEN
149 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
150 0 : matrix_type=dbcsr_type_no_symmetry)
151 : CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, &
152 0 : 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.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.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
159 :
160 0 : frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp2)
161 0 : CALL dbcsr_add_on_diag(matrix_tmp2, -1.0_dp)
162 0 : frob_matrix = dbcsr_frobenius_norm(matrix_tmp2)
163 0 : IF (unit_nr > 0) THEN
164 0 : WRITE (unit_nr, *) "Error for (inv(sqrt(S))*S*inv(sqrt(S))-I)", frob_matrix/frob_matrix_base
165 : END IF
166 :
167 0 : CALL dbcsr_release(matrix_tmp1)
168 0 : CALL dbcsr_release(matrix_tmp2)
169 : END IF
170 : END IF
171 :
172 : ! compute the inverse of S
173 12822 : 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 12792 : matrix_type=dbcsr_type_no_symmetry)
176 12792 : IF (.NOT. ls_scf_env%use_s_sqrt) THEN
177 2 : 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 12790 : 0.0_dp, ls_scf_env%matrix_s_inv, filter_eps=ls_scf_env%eps_filter)
181 : END IF
182 12792 : IF (ls_scf_env%check_s_inv) THEN
183 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
184 0 : 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.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
187 0 : frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp1)
188 0 : CALL dbcsr_add_on_diag(matrix_tmp1, -1.0_dp)
189 0 : frob_matrix = dbcsr_frobenius_norm(matrix_tmp1)
190 0 : IF (unit_nr > 0) THEN
191 0 : WRITE (unit_nr, *) "Error for (inv(S)*S-I)", frob_matrix/frob_matrix_base
192 : END IF
193 0 : CALL dbcsr_release(matrix_tmp1)
194 : END IF
195 : END IF
196 :
197 12822 : CALL timestop(handle)
198 12822 : 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 456 : 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 456 : REAL(dp), DIMENSION(:, :), POINTER :: block_dp
234 : TYPE(dbcsr_iterator_type) :: iter
235 : TYPE(dbcsr_type) :: matrix_bs
236 :
237 456 : CALL timeset(routineN, handle)
238 :
239 : ! first generate a block diagonal copy of s
240 456 : CALL dbcsr_create(matrix_bs, template=matrix_s)
241 :
242 912 : SELECT CASE (preconditioner_type)
243 : CASE (ls_s_preconditioner_none)
244 : CASE (ls_s_preconditioner_atomic, ls_s_preconditioner_molecular)
245 456 : CALL dbcsr_iterator_start(iter, matrix_s)
246 28479 : DO WHILE (dbcsr_iterator_blocks_left(iter))
247 28023 : 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 28023 : block_needed = .FALSE.
255 :
256 28023 : IF (iblock_row == iblock_col) THEN
257 : block_needed = .TRUE.
258 : ELSE
259 24925 : IF (preconditioner_type == ls_s_preconditioner_molecular .AND. &
260 : ls_mstruct%cluster_type == ls_cluster_atomic) THEN
261 5151 : 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 456 : IF (block_needed) THEN
267 4286 : CALL dbcsr_put_block(matrix=matrix_bs, row=iblock_row, col=iblock_col, block=block_dp)
268 : END IF
269 :
270 : END DO
271 912 : CALL dbcsr_iterator_stop(iter)
272 : END SELECT
273 :
274 456 : CALL dbcsr_finalize(matrix_bs)
275 :
276 456 : SELECT CASE (preconditioner_type)
277 : CASE (ls_s_preconditioner_none)
278 : ! for now make it a simple identity matrix
279 0 : CALL dbcsr_copy(matrix_bs_sqrt, matrix_bs)
280 0 : CALL dbcsr_set(matrix_bs_sqrt, 0.0_dp)
281 0 : CALL dbcsr_add_on_diag(matrix_bs_sqrt, 1.0_dp)
282 :
283 : ! for now make it a simple identity matrix
284 0 : CALL dbcsr_copy(matrix_bs_sqrt_inv, matrix_bs)
285 0 : CALL dbcsr_set(matrix_bs_sqrt_inv, 0.0_dp)
286 0 : CALL dbcsr_add_on_diag(matrix_bs_sqrt_inv, 1.0_dp)
287 : CASE (ls_s_preconditioner_atomic, ls_s_preconditioner_molecular)
288 456 : CALL dbcsr_copy(matrix_bs_sqrt, matrix_bs)
289 456 : 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 912 : iounit=-1)
298 : END SELECT
299 :
300 456 : CALL dbcsr_release(matrix_bs)
301 :
302 456 : CALL timestop(handle)
303 :
304 456 : 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 4804 : 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 4804 : CALL timeset(routineN, handle)
330 4804 : CALL dbcsr_create(matrix_tmp, template=matrix, matrix_type=dbcsr_type_no_symmetry)
331 :
332 4204 : SELECT CASE (direction)
333 : CASE ("forward")
334 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_bs_sqrt_inv, &
335 4204 : 0.0_dp, matrix_tmp)
336 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_bs_sqrt_inv, matrix_tmp, &
337 4204 : 0.0_dp, matrix)
338 : CASE ("backward")
339 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_bs_sqrt, &
340 600 : 0.0_dp, matrix_tmp)
341 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_bs_sqrt, matrix_tmp, &
342 600 : 0.0_dp, matrix)
343 : CASE DEFAULT
344 4804 : CPABORT("Direction should be forward or backward when applying preconditioner")
345 : END SELECT
346 :
347 4804 : CALL dbcsr_release(matrix_tmp)
348 :
349 4804 : CALL timestop(handle)
350 :
351 4804 : 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 1044 : 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 1044 : CALL timeset(routineN, handle)
400 :
401 1044 : logger => cp_get_default_logger()
402 1044 : IF (logger%para_env%is_source()) THEN
403 522 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
404 : ELSE
405 : unit_nr = -1
406 : END IF
407 :
408 1044 : do_sign_symmetric = .FALSE.
409 1044 : IF (PRESENT(sign_symmetric)) do_sign_symmetric = sign_symmetric
410 :
411 1044 : used_submatrix_sign_method = ls_scf_submatrix_sign_ns
412 1044 : 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 1044 : used_submatrix_sign_method == ls_scf_submatrix_sign_direct_muadj_lowmem))
417 :
418 4 : 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 4 : nelectron, matrix_s_sqrt_inv)
423 : ELSE
424 1040 : increment = initial_increment
425 :
426 1040 : has_mu_low = .FALSE.
427 1040 : has_mu_high = .FALSE.
428 :
429 : ! bisect if both bounds are known, otherwise find the bounds with a linear search
430 1184 : DO iter = 1, 30
431 1184 : IF (has_mu_low .AND. has_mu_high) THEN
432 16 : mu = (mu_low + mu_high)/2
433 16 : 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 1184 : matrix_s_sqrt_inv)
440 1184 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,I2,1X,F13.9,1X,F15.9)') &
441 592 : "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 1184 : IF (ABS(trace - nelectron) < 0.5_dp .OR. fixed_mu) EXIT
446 :
447 2368 : IF (trace < nelectron) THEN
448 32 : mu_low = mu
449 32 : mu = mu + increment
450 32 : has_mu_low = .TRUE.
451 32 : increment = increment*2
452 : ELSE
453 112 : mu_high = mu
454 112 : mu = mu - increment
455 112 : has_mu_high = .TRUE.
456 112 : increment = increment*2
457 : END IF
458 : END DO
459 :
460 : END IF
461 :
462 1044 : CALL timestop(handle)
463 :
464 1044 : 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 2412 : 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 1206 : CALL timeset(routineN, handle)
507 :
508 1206 : logger => cp_get_default_logger()
509 1206 : IF (logger%para_env%is_source()) THEN
510 603 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
511 : ELSE
512 : unit_nr = -1
513 : END IF
514 :
515 1206 : CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
516 :
517 1206 : IF (sign_symmetric) THEN
518 :
519 4 : IF (.NOT. PRESENT(matrix_s_sqrt_inv)) THEN
520 0 : CPABORT("Argument matrix_s_sqrt_inv required if sign_symmetric is set")
521 : END IF
522 :
523 4 : CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
524 4 : CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv2, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
525 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
526 4 : 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, filter_eps=threshold)
527 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, matrix_s_sqrt_inv, &
528 4 : 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
529 4 : CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)
530 :
531 6 : SELECT CASE (sign_method)
532 : CASE (ls_scf_sign_ns)
533 2 : CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, iounit=-1)
534 : CASE (ls_scf_sign_proot)
535 0 : CALL matrix_sign_proot(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
536 : CASE (ls_scf_sign_submatrix)
537 2 : CALL matrix_sign_submatrix(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, submatrix_sign_method)
538 : CASE DEFAULT
539 4 : CPABORT("Unkown sign method.")
540 : END SELECT
541 4 : CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
542 4 : CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2)
543 :
544 : ELSE ! .NOT. sign_symmetric
545 : ! get inv(S)*H-I*mu
546 1202 : CALL dbcsr_create(matrix_sinv_ks, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
547 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_inv, matrix_ks, &
548 1202 : 0.0_dp, matrix_sinv_ks, filter_eps=threshold)
549 1202 : CALL dbcsr_add_on_diag(matrix_sinv_ks, -mu)
550 :
551 : ! compute sign(inv(S)*H-I*mu)
552 2392 : SELECT CASE (sign_method)
553 : CASE (ls_scf_sign_ns)
554 1190 : CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order, iounit=-1)
555 : CASE (ls_scf_sign_proot)
556 8 : CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order)
557 : CASE (ls_scf_sign_submatrix)
558 4 : CALL matrix_sign_submatrix(matrix_sign, matrix_sinv_ks, threshold, sign_order, submatrix_sign_method)
559 : CASE DEFAULT
560 1202 : CPABORT("Unkown sign method.")
561 : END SELECT
562 1202 : CALL dbcsr_release(matrix_sinv_ks)
563 : END IF
564 :
565 : ! now construct the density matrix PS=0.5*(I-sign(inv(S)H-I*mu))
566 1206 : CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
567 1206 : CALL dbcsr_copy(matrix_p_ud, matrix_sign)
568 1206 : CALL dbcsr_scale(matrix_p_ud, -0.5_dp)
569 1206 : CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp)
570 1206 : CALL dbcsr_release(matrix_sign)
571 :
572 : ! we now have PS, lets get its trace
573 1206 : CALL dbcsr_trace(matrix_p_ud, trace)
574 :
575 : ! we can also check it is idempotent PS*PS=PS
576 1206 : CALL dbcsr_create(matrix_tmp, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
577 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_p_ud, &
578 1206 : 0.0_dp, matrix_tmp, filter_eps=threshold)
579 1206 : CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp)
580 1206 : frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
581 1206 : IF (unit_nr > 0 .AND. frob_matrix > 0.001_dp) THEN
582 30 : WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix
583 : END IF
584 :
585 1206 : IF (sign_symmetric) THEN
586 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, &
587 4 : 0.0_dp, matrix_tmp, filter_eps=threshold)
588 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_s_sqrt_inv, &
589 4 : 0.0_dp, matrix_p, filter_eps=threshold)
590 : ELSE
591 :
592 : ! get P=PS*inv(S)
593 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_s_inv, &
594 1202 : 0.0_dp, matrix_p, filter_eps=threshold)
595 : END IF
596 1206 : CALL dbcsr_release(matrix_p_ud)
597 1206 : CALL dbcsr_release(matrix_tmp)
598 :
599 1206 : CALL timestop(handle)
600 :
601 1206 : END SUBROUTINE density_matrix_sign_fixed_mu
602 :
603 : ! **************************************************************************************************
604 : !> \brief compute the corresponding density matrix and its trace, using methods with internal mu adjustment
605 : !> \param matrix_p ...
606 : !> \param trace ...
607 : !> \param mu ...
608 : !> \param sign_method ...
609 : !> \param matrix_ks ...
610 : !> \param matrix_s ...
611 : !> \param threshold ...
612 : !> \param submatrix_sign_method ...
613 : !> \param nelectron ...
614 : !> \param matrix_s_sqrt_inv ...
615 : !> \par History
616 : !> 2020.07 created, based on density_matrix_sign_fixed_mu [Michael Lass]
617 : !> \author Michael Lass
618 : ! **************************************************************************************************
619 8 : SUBROUTINE density_matrix_sign_internal_mu(matrix_p, trace, mu, sign_method, matrix_ks, &
620 : matrix_s, threshold, submatrix_sign_method, &
621 : nelectron, matrix_s_sqrt_inv)
622 :
623 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
624 : REAL(KIND=dp), INTENT(OUT) :: trace
625 : REAL(KIND=dp), INTENT(INOUT) :: mu
626 : INTEGER :: sign_method
627 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_ks, matrix_s
628 : REAL(KIND=dp), INTENT(IN) :: threshold
629 : INTEGER :: submatrix_sign_method
630 : INTEGER, INTENT(IN) :: nelectron
631 : TYPE(dbcsr_type), INTENT(IN) :: matrix_s_sqrt_inv
632 :
633 : CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_sign_internal_mu'
634 :
635 : INTEGER :: handle, unit_nr
636 : REAL(KIND=dp) :: frob_matrix
637 : TYPE(cp_logger_type), POINTER :: logger
638 : TYPE(dbcsr_type) :: matrix_p_ud, matrix_sign, &
639 : matrix_ssqrtinv_ks_ssqrtinv, &
640 : matrix_ssqrtinv_ks_ssqrtinv2, &
641 : matrix_tmp
642 :
643 4 : CALL timeset(routineN, handle)
644 :
645 4 : logger => cp_get_default_logger()
646 4 : IF (logger%para_env%is_source()) THEN
647 2 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
648 : ELSE
649 : unit_nr = -1
650 : END IF
651 :
652 4 : CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
653 :
654 4 : CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
655 4 : CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv2, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
656 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
657 4 : 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, filter_eps=threshold)
658 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_ssqrtinv_ks_ssqrtinv2, matrix_s_sqrt_inv, &
659 4 : 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
660 4 : CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)
661 :
662 8 : SELECT CASE (sign_method)
663 : CASE (ls_scf_sign_submatrix)
664 8 : SELECT CASE (submatrix_sign_method)
665 : CASE (ls_scf_submatrix_sign_direct_muadj, ls_scf_submatrix_sign_direct_muadj_lowmem)
666 : CALL matrix_sign_submatrix_mu_adjust(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, mu, nelectron, threshold, &
667 4 : submatrix_sign_method)
668 : CASE DEFAULT
669 4 : CPABORT("density_matrix_sign_internal_mu called with invalid submatrix sign method")
670 : END SELECT
671 : CASE DEFAULT
672 4 : CPABORT("density_matrix_sign_internal_mu called with invalid sign method.")
673 : END SELECT
674 4 : CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
675 4 : CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2)
676 :
677 : ! now construct the density matrix PS=0.5*(I-sign(inv(S)H-I*mu))
678 4 : CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
679 4 : CALL dbcsr_copy(matrix_p_ud, matrix_sign)
680 4 : CALL dbcsr_scale(matrix_p_ud, -0.5_dp)
681 4 : CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp)
682 4 : CALL dbcsr_release(matrix_sign)
683 :
684 : ! we now have PS, lets get its trace
685 4 : CALL dbcsr_trace(matrix_p_ud, trace)
686 :
687 : ! we can also check it is idempotent PS*PS=PS
688 4 : CALL dbcsr_create(matrix_tmp, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
689 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p_ud, matrix_p_ud, &
690 4 : 0.0_dp, matrix_tmp, filter_eps=threshold)
691 4 : CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp)
692 4 : frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
693 4 : IF (unit_nr > 0 .AND. frob_matrix > 0.001_dp) THEN
694 0 : WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix
695 : END IF
696 :
697 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, &
698 4 : 0.0_dp, matrix_tmp, filter_eps=threshold)
699 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_s_sqrt_inv, &
700 4 : 0.0_dp, matrix_p, filter_eps=threshold)
701 4 : CALL dbcsr_release(matrix_p_ud)
702 4 : CALL dbcsr_release(matrix_tmp)
703 :
704 4 : CALL timestop(handle)
705 :
706 4 : END SUBROUTINE density_matrix_sign_internal_mu
707 :
708 : ! **************************************************************************************************
709 : !> \brief compute the density matrix using a trace-resetting algorithm
710 : !> \param matrix_p ...
711 : !> \param matrix_ks ...
712 : !> \param matrix_s_sqrt_inv ...
713 : !> \param nelectron ...
714 : !> \param threshold ...
715 : !> \param e_homo ...
716 : !> \param e_lumo ...
717 : !> \param e_mu ...
718 : !> \param dynamic_threshold ...
719 : !> \param matrix_ks_deviation ...
720 : !> \param max_iter_lanczos ...
721 : !> \param eps_lanczos ...
722 : !> \param converged ...
723 : !> \param iounit ...
724 : !> \par History
725 : !> 2012.06 created [Florian Thoele]
726 : !> \author Florian Thoele
727 : ! **************************************************************************************************
728 13654 : SUBROUTINE density_matrix_trs4(matrix_p, matrix_ks, matrix_s_sqrt_inv, &
729 : nelectron, threshold, e_homo, e_lumo, e_mu, &
730 : dynamic_threshold, matrix_ks_deviation, &
731 : max_iter_lanczos, eps_lanczos, converged, iounit)
732 :
733 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
734 : TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv
735 : INTEGER, INTENT(IN) :: nelectron
736 : REAL(KIND=dp), INTENT(IN) :: threshold
737 : REAL(KIND=dp), INTENT(INOUT) :: e_homo, e_lumo, e_mu
738 : LOGICAL, INTENT(IN), OPTIONAL :: dynamic_threshold
739 : TYPE(dbcsr_type), INTENT(INOUT), OPTIONAL :: matrix_ks_deviation
740 : INTEGER, INTENT(IN) :: max_iter_lanczos
741 : REAL(KIND=dp), INTENT(IN) :: eps_lanczos
742 : LOGICAL, INTENT(OUT), OPTIONAL :: converged
743 : INTEGER, INTENT(IN), OPTIONAL :: iounit
744 :
745 : CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_trs4'
746 : INTEGER, PARAMETER :: max_iter = 100
747 : REAL(KIND=dp), PARAMETER :: gamma_max = 6.0_dp, gamma_min = 0.0_dp
748 :
749 : INTEGER :: branch, estimated_steps, handle, i, j, &
750 : unit_nr
751 : INTEGER(kind=int_8) :: flop1, flop2
752 : LOGICAL :: arnoldi_converged, do_dyn_threshold
753 : REAL(KIND=dp) :: current_threshold, delta_n, eps_max, eps_min, est_threshold, frob_id, &
754 : frob_x, gam, homo, lumo, max_eig, max_threshold, maxdev, maxev, min_eig, minev, mmin, mu, &
755 : mu_a, mu_b, mu_c, mu_fa, mu_fc, occ_matrix, scaled_homo_bound, scaled_lumo_bound, t1, t2, &
756 : trace_fx, trace_gx, xi
757 13654 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: gamma_values
758 : TYPE(cp_logger_type), POINTER :: logger
759 : TYPE(dbcsr_type) :: matrix_k0, matrix_x, matrix_x_nosym, &
760 : matrix_xidsq, matrix_xsq, tmp_gx
761 :
762 13654 : IF (nelectron == 0) THEN
763 0 : CALL dbcsr_set(matrix_p, 0.0_dp)
764 : RETURN
765 : END IF
766 :
767 13654 : CALL timeset(routineN, handle)
768 :
769 13654 : IF (PRESENT(iounit)) THEN
770 1818 : unit_nr = iounit
771 : ELSE
772 11836 : logger => cp_get_default_logger()
773 11836 : IF (logger%para_env%is_source()) THEN
774 5918 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
775 : ELSE
776 5918 : unit_nr = -1
777 : END IF
778 : END IF
779 :
780 13654 : do_dyn_threshold = .FALSE.
781 13654 : IF (PRESENT(dynamic_threshold)) do_dyn_threshold = dynamic_threshold
782 :
783 13654 : IF (PRESENT(converged)) converged = .FALSE.
784 :
785 : ! init X = (eps_n*I - H)/(eps_n - eps_0) ... H* = S^-1/2*H*S^-1/2
786 13654 : CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type="S")
787 :
788 : ! at some points the non-symmetric version of x is required
789 13654 : CALL dbcsr_create(matrix_x_nosym, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
790 :
791 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
792 13654 : 0.0_dp, matrix_x_nosym, filter_eps=threshold)
793 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x_nosym, matrix_s_sqrt_inv, &
794 13654 : 0.0_dp, matrix_x, filter_eps=threshold)
795 13654 : CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym)
796 :
797 13654 : CALL dbcsr_create(matrix_k0, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
798 13654 : CALL dbcsr_copy(matrix_k0, matrix_x_nosym)
799 :
800 : ! compute the deviation in the mixed matrix, as seen in the ortho basis
801 13654 : IF (do_dyn_threshold) THEN
802 24 : CPASSERT(PRESENT(matrix_ks_deviation))
803 24 : CALL dbcsr_add(matrix_ks_deviation, matrix_x_nosym, -1.0_dp, 1.0_dp)
804 : CALL arnoldi_extremal(matrix_ks_deviation, maxev, minev, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
805 24 : converged=arnoldi_converged)
806 24 : maxdev = MAX(ABS(maxev), ABS(minev))
807 24 : IF (unit_nr > 0) THEN
808 0 : WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", arnoldi_converged
809 0 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "change in mixed matrix: ", maxdev
810 0 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo + maxdev
811 0 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo - maxdev
812 0 : WRITE (unit_nr, '(T6,A,1X,L12)') "Predicts a gap ? ", ((e_lumo - maxdev) - (e_homo + maxdev)) > 0
813 : END IF
814 : ! save the old mixed matrix
815 24 : CALL dbcsr_copy(matrix_ks_deviation, matrix_x_nosym)
816 :
817 : END IF
818 :
819 : ! get largest/smallest eigenvalues for scaling
820 : CALL arnoldi_extremal(matrix_x_nosym, max_eig, min_eig, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
821 13654 : converged=arnoldi_converged)
822 19572 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", &
823 11836 : min_eig, max_eig, " converged: ", arnoldi_converged
824 13654 : eps_max = max_eig
825 13654 : eps_min = min_eig
826 :
827 : ! scale KS matrix
828 13654 : IF (eps_max == eps_min) THEN
829 20 : CALL dbcsr_scale(matrix_x, 1.0_dp/eps_max)
830 : ELSE
831 13634 : CALL dbcsr_add_on_diag(matrix_x, -eps_max)
832 13634 : CALL dbcsr_scale(matrix_x, -1.0_dp/(eps_max - eps_min))
833 : END IF
834 :
835 13654 : current_threshold = threshold
836 13654 : IF (do_dyn_threshold) THEN
837 : ! scale bounds for HOMO/LUMO
838 24 : scaled_homo_bound = (eps_max - (e_homo + maxdev))/(eps_max - eps_min)
839 24 : scaled_lumo_bound = (eps_max - (e_lumo - maxdev))/(eps_max - eps_min)
840 : END IF
841 :
842 13654 : CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type="S")
843 :
844 13654 : CALL dbcsr_create(matrix_xidsq, template=matrix_ks, matrix_type="S")
845 :
846 13654 : CALL dbcsr_create(tmp_gx, template=matrix_ks, matrix_type="S")
847 :
848 13654 : ALLOCATE (gamma_values(max_iter))
849 :
850 69948 : DO i = 1, max_iter
851 69948 : t1 = m_walltime()
852 69948 : flop1 = 0; flop2 = 0
853 :
854 : ! get X*X
855 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, &
856 : 0.0_dp, matrix_xsq, &
857 69948 : filter_eps=current_threshold, flop=flop1)
858 :
859 : ! intermediate use matrix_xidsq to compute = X*X-X
860 69948 : CALL dbcsr_copy(matrix_xidsq, matrix_x)
861 69948 : CALL dbcsr_add(matrix_xidsq, matrix_xsq, -1.0_dp, 1.0_dp)
862 69948 : frob_id = dbcsr_frobenius_norm(matrix_xidsq)
863 69948 : frob_x = dbcsr_frobenius_norm(matrix_x)
864 :
865 : ! xidsq = (1-X)*(1-X)
866 : ! use (1-x)*(1-x) = 1 + x*x - 2*x
867 69948 : CALL dbcsr_copy(matrix_xidsq, matrix_x)
868 69948 : CALL dbcsr_add(matrix_xidsq, matrix_xsq, -2.0_dp, 1.0_dp)
869 69948 : CALL dbcsr_add_on_diag(matrix_xidsq, 1.0_dp)
870 :
871 : ! tmp_gx = 4X-3X*X
872 69948 : CALL dbcsr_copy(tmp_gx, matrix_x)
873 69948 : CALL dbcsr_add(tmp_gx, matrix_xsq, 4.0_dp, -3.0_dp)
874 :
875 : ! get gamma
876 : ! Tr(F) = Tr(XX*tmp_gx) Tr(G) is equivalent
877 69948 : CALL dbcsr_dot(matrix_xsq, matrix_xidsq, trace_gx)
878 69948 : CALL dbcsr_dot(matrix_xsq, tmp_gx, trace_fx)
879 :
880 : ! if converged, and gam becomes noisy, fix it to 3, which results in a final McWeeny step.
881 : ! do this only if the electron count is reasonable.
882 : ! maybe tune if the current criterion is not good enough
883 69948 : delta_n = nelectron - trace_fx
884 : ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ...
885 69948 : IF (((frob_id*frob_id) < (threshold*frob_x*frob_x)) .AND. (ABS(delta_n) < 0.5_dp)) THEN
886 13654 : gam = 3.0_dp
887 56294 : ELSE IF (ABS(delta_n) < 1e-14_dp) THEN
888 0 : gam = 0.0_dp ! rare case of perfect electron count
889 : ELSE
890 : ! make sure, we don't divide by zero, as soon as gam is outside the interval gam_min,gam_max, it doesn't matter
891 56294 : gam = delta_n/MAX(trace_gx, ABS(delta_n)/100)
892 : END IF
893 69948 : gamma_values(i) = gam
894 :
895 : IF (unit_nr > 0 .AND. .FALSE.) THEN
896 : WRITE (unit_nr, *) "trace_fx", trace_fx, "trace_gx", trace_gx, "gam", gam, &
897 : "frob_id", frob_id, "conv", ABS(frob_id/frob_x)
898 : END IF
899 :
900 69948 : IF (do_dyn_threshold) THEN
901 : ! quantities used for dynamic thresholding, when the estimated gap is larger than zero
902 154 : xi = (scaled_homo_bound - scaled_lumo_bound)
903 154 : IF (xi > 0.0_dp) THEN
904 130 : mmin = 0.5*(scaled_homo_bound + scaled_lumo_bound)
905 130 : max_threshold = ABS(1 - 2*mmin)*xi
906 :
907 130 : scaled_homo_bound = evaluate_trs4_polynomial(scaled_homo_bound, gamma_values(i:), 1)
908 130 : scaled_lumo_bound = evaluate_trs4_polynomial(scaled_lumo_bound, gamma_values(i:), 1)
909 130 : estimated_steps = estimate_steps(scaled_homo_bound, scaled_lumo_bound, threshold)
910 :
911 130 : est_threshold = (threshold/(estimated_steps + i + 1))*xi/(1 + threshold/(estimated_steps + i + 1))
912 130 : est_threshold = MIN(max_threshold, est_threshold)
913 130 : IF (i > 1) est_threshold = MAX(est_threshold, 0.1_dp*current_threshold)
914 130 : current_threshold = est_threshold
915 : ELSE
916 24 : current_threshold = threshold
917 : END IF
918 : END IF
919 :
920 69948 : IF (gam > gamma_max) THEN
921 : ! Xn+1 = 2X-X*X
922 700 : CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp)
923 700 : CALL dbcsr_filter(matrix_x, current_threshold)
924 700 : branch = 1
925 69248 : ELSE IF (gam < gamma_min) THEN
926 : ! Xn+1 = X*X
927 3594 : CALL dbcsr_copy(matrix_x, matrix_xsq)
928 3594 : branch = 2
929 : ELSE
930 : ! Xn+1 = F(X) + gam*G(X)
931 65654 : CALL dbcsr_add(tmp_gx, matrix_xidsq, 1.0_dp, gam)
932 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_xsq, tmp_gx, &
933 : 0.0_dp, matrix_x, &
934 65654 : flop=flop2, filter_eps=current_threshold)
935 65654 : branch = 3
936 : END IF
937 :
938 69948 : occ_matrix = dbcsr_get_occupation(matrix_x)
939 69948 : t2 = m_walltime()
940 69948 : IF (unit_nr > 0) THEN
941 : WRITE (unit_nr, &
942 30176 : '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TRS4 it ", &
943 30176 : i, occ_matrix, ABS(trace_gx), t2 - t1, &
944 60352 : (flop1 + flop2)/(1.0E6_dp*MAX(t2 - t1, 0.001_dp)), current_threshold
945 30176 : CALL m_flush(unit_nr)
946 : END IF
947 :
948 69948 : IF (abnormal_value(trace_gx)) THEN
949 0 : CPABORT("trace_gx is an abnormal value (NaN/Inf).")
950 : END IF
951 :
952 : ! a branch of 1 or 2 appears to lead to a less accurate electron number count and premature exit
953 : ! if it turns out this does not exit because we get stuck in branch 1/2 for a reason we need to refine further
954 : ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ...
955 139896 : IF ((frob_id*frob_id) < (threshold*frob_x*frob_x) .AND. branch == 3 .AND. (ABS(delta_n) < 0.5_dp)) THEN
956 13654 : IF (PRESENT(converged)) converged = .TRUE.
957 : EXIT
958 : END IF
959 :
960 : END DO
961 :
962 13654 : occ_matrix = dbcsr_get_occupation(matrix_x)
963 13654 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,I3,1X,F10.8,E12.3)') 'Final TRS4 iteration ', i, occ_matrix, ABS(trace_gx)
964 :
965 : ! free some memory
966 13654 : CALL dbcsr_release(tmp_gx)
967 13654 : CALL dbcsr_release(matrix_xsq)
968 13654 : CALL dbcsr_release(matrix_xidsq)
969 :
970 : ! output to matrix_p, P = inv(S)^0.5 X inv(S)^0.5
971 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, &
972 13654 : 0.0_dp, matrix_x_nosym, filter_eps=threshold)
973 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_x_nosym, &
974 13654 : 0.0_dp, matrix_p, filter_eps=threshold)
975 :
976 : ! calculate the chemical potential by doing a bisection of fk(x0)-0.5,
977 : ! where fk is evaluated using the stored values for gamma
978 : ! E. Rubensson et al., Chem Phys Lett 432, 2006, 591-594
979 13654 : mu_a = 0.0_dp; mu_b = 1.0_dp
980 13654 : mu_fa = evaluate_trs4_polynomial(mu_a, gamma_values, i - 1) - 0.5_dp
981 160868 : DO j = 1, 40
982 160868 : mu_c = 0.5*(mu_a + mu_b)
983 : ! i-1 because in the last iteration, only convergence is checked
984 160868 : mu_fc = evaluate_trs4_polynomial(mu_c, gamma_values, i - 1) - 0.5_dp
985 160868 : IF (ABS(mu_fc) < 1.0E-6_dp .OR. (mu_b - mu_a)/2 < 1.0E-6_dp) EXIT !TODO: define threshold values
986 :
987 160868 : IF (mu_fc*mu_fa > 0) THEN
988 75942 : mu_a = mu_c
989 75942 : mu_fa = mu_fc
990 : ELSE
991 : mu_b = mu_c
992 : END IF
993 : END DO
994 13654 : mu = (eps_min - eps_max)*mu_c + eps_max
995 13654 : DEALLOCATE (gamma_values)
996 13654 : IF (unit_nr > 0) THEN
997 5918 : WRITE (unit_nr, '(T6,A,1X,F12.5)') 'Chemical potential (mu): ', mu
998 : END IF
999 13654 : e_mu = mu
1000 :
1001 13654 : IF (do_dyn_threshold) THEN
1002 24 : CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym)
1003 : CALL compute_homo_lumo(matrix_k0, matrix_x_nosym, eps_min, eps_max, &
1004 24 : threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr)
1005 24 : e_homo = homo
1006 24 : e_lumo = lumo
1007 : END IF
1008 :
1009 13654 : CALL dbcsr_release(matrix_x)
1010 13654 : CALL dbcsr_release(matrix_x_nosym)
1011 13654 : CALL dbcsr_release(matrix_k0)
1012 13654 : CALL timestop(handle)
1013 :
1014 27308 : END SUBROUTINE density_matrix_trs4
1015 :
1016 : ! **************************************************************************************************
1017 : !> \brief compute the density matrix using a non monotonic trace conserving
1018 : !> algorithm based on SIAM DOI. 10.1137/130911585.
1019 : !> 2014.04 created [Jonathan Mullin]
1020 : !> \param matrix_p ...
1021 : !> \param matrix_ks ...
1022 : !> \param matrix_s_sqrt_inv ...
1023 : !> \param nelectron ...
1024 : !> \param threshold ...
1025 : !> \param e_homo ...
1026 : !> \param e_lumo ...
1027 : !> \param non_monotonic ...
1028 : !> \param eps_lanczos ...
1029 : !> \param max_iter_lanczos ...
1030 : !> \param iounit ...
1031 : !> \author Jonathan Mullin
1032 : ! **************************************************************************************************
1033 286 : SUBROUTINE density_matrix_tc2(matrix_p, matrix_ks, matrix_s_sqrt_inv, &
1034 : nelectron, threshold, e_homo, e_lumo, &
1035 : non_monotonic, eps_lanczos, max_iter_lanczos, iounit)
1036 :
1037 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_p
1038 : TYPE(dbcsr_type), INTENT(IN) :: matrix_ks, matrix_s_sqrt_inv
1039 : INTEGER, INTENT(IN) :: nelectron
1040 : REAL(KIND=dp), INTENT(IN) :: threshold
1041 : REAL(KIND=dp), INTENT(INOUT) :: e_homo, e_lumo
1042 : LOGICAL, INTENT(IN), OPTIONAL :: non_monotonic
1043 : REAL(KIND=dp), INTENT(IN) :: eps_lanczos
1044 : INTEGER, INTENT(IN) :: max_iter_lanczos
1045 : INTEGER, INTENT(IN), OPTIONAL :: iounit
1046 :
1047 : CHARACTER(LEN=*), PARAMETER :: routineN = 'density_matrix_tc2'
1048 : INTEGER, PARAMETER :: max_iter = 100
1049 :
1050 : INTEGER :: handle, i, j, k, unit_nr
1051 : INTEGER(kind=int_8) :: flop1, flop2
1052 : LOGICAL :: converged, do_non_monotonic, &
1053 : tc2_converged
1054 : REAL(KIND=dp) :: beta, betaB, eps_max, eps_min, gama, &
1055 : max_eig, min_eig, occ_matrix, t1, t2, &
1056 : trace_fx, trace_gx
1057 286 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, lambda, nu, poly, wu, X, Y
1058 : TYPE(cp_logger_type), POINTER :: logger
1059 : TYPE(dbcsr_type) :: matrix_tmp, matrix_x, matrix_xsq
1060 :
1061 286 : CALL timeset(routineN, handle)
1062 :
1063 286 : IF (PRESENT(iounit)) THEN
1064 284 : unit_nr = iounit
1065 : ELSE
1066 2 : logger => cp_get_default_logger()
1067 2 : IF (logger%para_env%is_source()) THEN
1068 1 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1069 : ELSE
1070 1 : unit_nr = -1
1071 : END IF
1072 : END IF
1073 :
1074 286 : do_non_monotonic = .FALSE.
1075 286 : IF (PRESENT(non_monotonic)) do_non_monotonic = non_monotonic
1076 :
1077 : ! init X = (eps_n*I - H)/(eps_n - eps_0) ... H* = S^-1/2*H*S^-1/2
1078 286 : CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1079 286 : CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1080 :
1081 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_ks, &
1082 286 : 0.0_dp, matrix_xsq, filter_eps=threshold)
1083 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_xsq, matrix_s_sqrt_inv, &
1084 286 : 0.0_dp, matrix_x, filter_eps=threshold)
1085 :
1086 286 : IF (unit_nr > 0) THEN
1087 1 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo
1088 1 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo
1089 1 : WRITE (unit_nr, '(T6,A,1X,L12)') "Predicts a gap ? ", ((e_lumo) - (e_homo)) > 0
1090 : END IF
1091 :
1092 : ! get largest/smallest eigenvalues for scaling
1093 : CALL arnoldi_extremal(matrix_x, max_eig, min_eig, max_iter=max_iter_lanczos, threshold=eps_lanczos, &
1094 286 : converged=converged)
1095 287 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", &
1096 2 : min_eig, max_eig, " converged: ", converged
1097 :
1098 286 : eps_max = max_eig
1099 286 : eps_min = min_eig
1100 :
1101 : ! scale KS matrix
1102 286 : CALL dbcsr_scale(matrix_x, -1.0_dp)
1103 286 : CALL dbcsr_add_on_diag(matrix_x, eps_max)
1104 286 : CALL dbcsr_scale(matrix_x, 1/(eps_max - eps_min))
1105 :
1106 286 : CALL dbcsr_copy(matrix_xsq, matrix_x)
1107 :
1108 286 : CALL dbcsr_create(matrix_tmp, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1109 :
1110 286 : ALLOCATE (poly(max_iter))
1111 286 : ALLOCATE (nu(max_iter))
1112 286 : ALLOCATE (wu(max_iter))
1113 286 : ALLOCATE (alpha(max_iter))
1114 286 : ALLOCATE (X(4))
1115 286 : ALLOCATE (Y(4))
1116 286 : ALLOCATE (lambda(4))
1117 :
1118 : ! Controls over the non_monotonic bounds, First if low gap, bias slightly
1119 286 : beta = (eps_max - ABS(e_lumo))/(eps_max - eps_min)
1120 286 : betaB = (eps_max + ABS(e_homo))/(eps_max - eps_min)
1121 :
1122 286 : IF ((beta - betaB) < 0.005_dp) THEN
1123 286 : beta = beta - 0.002_dp
1124 286 : betaB = betaB + 0.002_dp
1125 : END IF
1126 : ! Check if input specifies to use monotonic bounds.
1127 286 : IF (.NOT. do_non_monotonic) THEN
1128 26 : beta = 0.0_dp
1129 26 : betaB = 1.0_dp
1130 : END IF
1131 : ! initial SCF cycle has no reliable estimate of homo/lumo, force monotinic bounds.
1132 286 : IF (e_homo == 0.0_dp) THEN
1133 102 : beta = 0.0_dp
1134 102 : BetaB = 1.0_dp
1135 : END IF
1136 :
1137 : ! init to take true branch first
1138 286 : trace_fx = nelectron
1139 286 : trace_gx = 0
1140 :
1141 286 : tc2_converged = .FALSE.
1142 4636 : DO i = 1, max_iter
1143 4636 : t1 = m_walltime()
1144 4636 : flop1 = 0; flop2 = 0
1145 :
1146 4636 : IF (ABS(trace_fx - nelectron) <= ABS(trace_gx - nelectron)) THEN
1147 : ! Xn+1 = (aX+ (1-a)I)^2
1148 2418 : poly(i) = 1.0_dp
1149 2418 : alpha(i) = 2.0_dp/(2.0_dp - beta)
1150 :
1151 2418 : CALL dbcsr_scale(matrix_x, alpha(i))
1152 2418 : CALL dbcsr_add_on_diag(matrix_x, 1.0_dp - alpha(i))
1153 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, &
1154 : 0.0_dp, matrix_xsq, &
1155 2418 : filter_eps=threshold, flop=flop1)
1156 :
1157 : !save X for control variables
1158 2418 : CALL dbcsr_copy(matrix_tmp, matrix_x)
1159 :
1160 2418 : CALL dbcsr_copy(matrix_x, matrix_xsq)
1161 :
1162 2418 : beta = (1.0_dp - alpha(i)) + alpha(i)*beta
1163 2418 : beta = beta*beta
1164 2418 : betaB = (1.0_dp - alpha(i)) + alpha(i)*betaB
1165 2418 : betaB = betaB*betaB
1166 : ELSE
1167 : ! Xn+1 = 2aX-a^2*X^2
1168 2218 : poly(i) = 0.0_dp
1169 2218 : alpha(i) = 2.0_dp/(1.0_dp + betaB)
1170 :
1171 2218 : CALL dbcsr_scale(matrix_x, alpha(i))
1172 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_x, &
1173 : 0.0_dp, matrix_xsq, &
1174 2218 : filter_eps=threshold, flop=flop1)
1175 :
1176 : !save X for control variables
1177 2218 : CALL dbcsr_copy(matrix_tmp, matrix_x)
1178 : !
1179 2218 : CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp)
1180 :
1181 2218 : beta = alpha(i)*beta
1182 2218 : beta = 2.0_dp*beta - beta*beta
1183 2218 : betaB = alpha(i)*betaB
1184 2218 : betaB = 2.0_dp*betaB - betaB*betaB
1185 :
1186 : END IF
1187 4636 : occ_matrix = dbcsr_get_occupation(matrix_x)
1188 4636 : t2 = m_walltime()
1189 4636 : IF (unit_nr > 0) THEN
1190 : WRITE (unit_nr, &
1191 18 : '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TC2 it ", &
1192 18 : i, occ_matrix, t2 - t1, &
1193 36 : (flop1 + flop2)/(1.0E6_dp*(t2 - t1)), threshold
1194 18 : CALL m_flush(unit_nr)
1195 : END IF
1196 :
1197 : ! calculate control terms
1198 4636 : CALL dbcsr_trace(matrix_xsq, trace_fx)
1199 :
1200 : ! intermediate use matrix_xsq compute X- X*X , temorarily use trace_gx
1201 4636 : CALL dbcsr_add(matrix_xsq, matrix_tmp, -1.0_dp, 1.0_dp)
1202 4636 : CALL dbcsr_trace(matrix_xsq, trace_gx)
1203 4636 : nu(i) = dbcsr_frobenius_norm(matrix_xsq)
1204 4636 : wu(i) = trace_gx
1205 :
1206 : ! intermediate use matrix_xsq to compute = 2X - X*X
1207 4636 : CALL dbcsr_add(matrix_xsq, matrix_tmp, 1.0_dp, 1.0_dp)
1208 4636 : CALL dbcsr_trace(matrix_xsq, trace_gx)
1209 : ! TC2 has quadratic convergence, using the frobeniums norm as an idempotency deviation test.
1210 13908 : IF (ABS(nu(i)) < (threshold)) THEN
1211 : tc2_converged = .TRUE.
1212 : EXIT
1213 : END IF
1214 : END DO
1215 286 : IF (.NOT. tc2_converged) i = max_iter
1216 :
1217 286 : occ_matrix = dbcsr_get_occupation(matrix_x)
1218 286 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,I3,1X,1F10.8,1X,1F10.8)') 'Final TC2 iteration ', i, occ_matrix, ABS(nu(i))
1219 :
1220 : ! output to matrix_p, P = inv(S)^0.5 X inv(S)^0.5
1221 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_x, matrix_s_sqrt_inv, &
1222 286 : 0.0_dp, matrix_tmp, filter_eps=threshold)
1223 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_tmp, &
1224 286 : 0.0_dp, matrix_p, filter_eps=threshold)
1225 :
1226 286 : CALL dbcsr_release(matrix_xsq)
1227 286 : CALL dbcsr_release(matrix_tmp)
1228 :
1229 : ! ALGO 3 from. SIAM DOI. 10.1137/130911585
1230 286 : X(1) = 1.0_dp
1231 286 : X(2) = 1.0_dp
1232 286 : X(3) = 0.0_dp
1233 286 : X(4) = 0.0_dp
1234 : gama = 6.0_dp - 4.0_dp*(SQRT(2.0_dp))
1235 286 : gama = gama - gama*gama
1236 : DO
1237 2698 : IF (i < 1) EXIT
1238 2698 : IF (nu(i) >= gama) EXIT
1239 : ! safeguard against negative root, is skipping correct?
1240 2412 : IF (wu(i) < 1.0e-14_dp) THEN
1241 42 : i = i - 1
1242 42 : CYCLE
1243 : END IF
1244 2370 : IF ((1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)) < 0.0_dp) THEN
1245 12 : i = i - 1
1246 12 : CYCLE
1247 : END IF
1248 2358 : Y(1) = 0.5_dp*(1.0_dp - SQRT(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)))
1249 2358 : Y(2) = 0.5_dp*(1.0_dp - SQRT(1.0_dp - 4.0_dp*nu(i)))
1250 2358 : Y(3) = 0.5_dp*(1.0_dp + SQRT(1.0_dp - 4.0_dp*nu(i)))
1251 2358 : Y(4) = 0.5_dp*(1.0_dp + SQRT(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)))
1252 11790 : Y(:) = MIN(1.0_dp, MAX(0.0_dp, Y(:)))
1253 33026 : DO j = i, 1, -1
1254 33026 : IF (poly(j) == 1.0_dp) THEN
1255 81120 : DO k = 1, 4
1256 64896 : Y(k) = SQRT(Y(k))
1257 81120 : Y(k) = (Y(k) - 1.0_dp + alpha(j))/alpha(j)
1258 : END DO ! end K
1259 : ELSE
1260 72220 : DO k = 1, 4
1261 57776 : Y(k) = 1.0_dp - SQRT(1.0_dp - Y(k))
1262 72220 : Y(k) = Y(k)/alpha(j)
1263 : END DO ! end K
1264 : END IF ! end poly
1265 : END DO ! end j
1266 2358 : X(1) = MIN(X(1), Y(1))
1267 2358 : X(2) = MIN(X(2), Y(2))
1268 2358 : X(3) = MAX(X(3), Y(3))
1269 2358 : X(4) = MAX(X(4), Y(4))
1270 2644 : i = i - 1
1271 : END DO ! end i
1272 : ! lambda 1,2,3,4 are:: out lumo, in lumo, in homo, out homo
1273 1430 : DO k = 1, 4
1274 1430 : lambda(k) = eps_max - (eps_max - eps_min)*X(k)
1275 : END DO ! end k
1276 : ! END ALGO 3 from. SIAM DOI. 10.1137/130911585
1277 286 : e_homo = lambda(4)
1278 286 : e_lumo = lambda(1)
1279 286 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "outer homo/lumo/gap", e_homo, e_lumo, (e_lumo - e_homo)
1280 286 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "inner homo/lumo/gap", lambda(3), lambda(2), (lambda(2) - lambda(3))
1281 :
1282 286 : DEALLOCATE (poly)
1283 286 : DEALLOCATE (nu)
1284 286 : DEALLOCATE (wu)
1285 286 : DEALLOCATE (alpha)
1286 286 : DEALLOCATE (X)
1287 286 : DEALLOCATE (Y)
1288 286 : DEALLOCATE (lambda)
1289 :
1290 286 : CALL dbcsr_release(matrix_x)
1291 286 : CALL timestop(handle)
1292 :
1293 572 : END SUBROUTINE density_matrix_tc2
1294 :
1295 : ! **************************************************************************************************
1296 : !> \brief compute the homo and lumo given a KS matrix and a density matrix in the orthonormalized basis
1297 : !> and the eps_min and eps_max, min and max eigenvalue of the ks matrix
1298 : !> \param matrix_k ...
1299 : !> \param matrix_p ...
1300 : !> \param eps_min ...
1301 : !> \param eps_max ...
1302 : !> \param threshold ...
1303 : !> \param max_iter_lanczos ...
1304 : !> \param eps_lanczos ...
1305 : !> \param homo ...
1306 : !> \param lumo ...
1307 : !> \param unit_nr ...
1308 : !> \par History
1309 : !> 2012.06 created [Florian Thoele]
1310 : !> \author Florian Thoele
1311 : ! **************************************************************************************************
1312 132 : SUBROUTINE compute_homo_lumo(matrix_k, matrix_p, eps_min, eps_max, threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr)
1313 : TYPE(dbcsr_type) :: matrix_k, matrix_p
1314 : REAL(KIND=dp) :: eps_min, eps_max, threshold
1315 : INTEGER, INTENT(IN) :: max_iter_lanczos
1316 : REAL(KIND=dp), INTENT(IN) :: eps_lanczos
1317 : REAL(KIND=dp) :: homo, lumo
1318 : INTEGER :: unit_nr
1319 :
1320 : LOGICAL :: converged
1321 : REAL(KIND=dp) :: max_eig, min_eig, shift1, shift2
1322 : TYPE(dbcsr_type) :: tmp1, tmp2, tmp3
1323 :
1324 : ! temporary matrices used for HOMO/LUMO calculation
1325 :
1326 44 : CALL dbcsr_create(tmp1, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1327 :
1328 44 : CALL dbcsr_create(tmp2, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1329 :
1330 44 : CALL dbcsr_create(tmp3, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1331 :
1332 44 : shift1 = -eps_min
1333 44 : shift2 = eps_max
1334 :
1335 : ! find largest ev of P*(K+shift*1), where shift is the neg. val. of the smallest ev of K
1336 44 : CALL dbcsr_copy(tmp2, matrix_k)
1337 44 : CALL dbcsr_add_on_diag(tmp2, shift1)
1338 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, tmp2, &
1339 44 : 0.0_dp, tmp1, filter_eps=threshold)
1340 : CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, &
1341 44 : threshold=eps_lanczos, max_iter=max_iter_lanczos)
1342 44 : homo = max_eig - shift1
1343 44 : IF (unit_nr > 0) THEN
1344 10 : WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged
1345 : END IF
1346 :
1347 : ! -(1-P)*(K-shift*1) = (1-P)*(shift*1 - K), where shift is the largest ev of K
1348 44 : CALL dbcsr_copy(tmp3, matrix_p)
1349 44 : CALL dbcsr_scale(tmp3, -1.0_dp)
1350 44 : CALL dbcsr_add_on_diag(tmp3, 1.0_dp) !tmp3 = 1-P
1351 44 : CALL dbcsr_copy(tmp2, matrix_k)
1352 44 : CALL dbcsr_add_on_diag(tmp2, -shift2)
1353 : CALL dbcsr_multiply("N", "N", -1.0_dp, tmp3, tmp2, &
1354 44 : 0.0_dp, tmp1, filter_eps=threshold)
1355 : CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, &
1356 44 : threshold=eps_lanczos, max_iter=max_iter_lanczos)
1357 44 : lumo = -max_eig + shift2
1358 :
1359 44 : IF (unit_nr > 0) THEN
1360 10 : WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged
1361 10 : WRITE (unit_nr, '(T6,A,1X,3F12.5)') 'HOMO/LUMO/gap', homo, lumo, lumo - homo
1362 : END IF
1363 44 : CALL dbcsr_release(tmp1)
1364 44 : CALL dbcsr_release(tmp2)
1365 44 : CALL dbcsr_release(tmp3)
1366 :
1367 44 : END SUBROUTINE compute_homo_lumo
1368 :
1369 : ! **************************************************************************************************
1370 : !> \brief ...
1371 : !> \param x ...
1372 : !> \param gamma_values ...
1373 : !> \param i ...
1374 : !> \return ...
1375 : ! **************************************************************************************************
1376 174782 : FUNCTION evaluate_trs4_polynomial(x, gamma_values, i) RESULT(xr)
1377 : REAL(KIND=dp) :: x
1378 : REAL(KIND=dp), DIMENSION(:) :: gamma_values
1379 : INTEGER :: i
1380 : REAL(KIND=dp) :: xr
1381 :
1382 : REAL(KIND=dp), PARAMETER :: gam_max = 6.0_dp, gam_min = 0.0_dp
1383 :
1384 : INTEGER :: k
1385 :
1386 174782 : xr = x
1387 1352062 : DO k = 1, i
1388 1352062 : IF (gamma_values(k) > gam_max) THEN
1389 14656 : xr = 2*xr - xr**2
1390 1162624 : ELSE IF (gamma_values(k) < gam_min) THEN
1391 75304 : xr = xr**2
1392 : ELSE
1393 1087320 : xr = (xr*xr)*(4*xr - 3*xr*xr) + gamma_values(k)*xr*xr*((1 - xr)**2)
1394 : END IF
1395 : END DO
1396 174782 : END FUNCTION evaluate_trs4_polynomial
1397 :
1398 : ! **************************************************************************************************
1399 : !> \brief ...
1400 : !> \param homo ...
1401 : !> \param lumo ...
1402 : !> \param threshold ...
1403 : !> \return ...
1404 : ! **************************************************************************************************
1405 130 : FUNCTION estimate_steps(homo, lumo, threshold) RESULT(steps)
1406 : REAL(KIND=dp) :: homo, lumo, threshold
1407 : INTEGER :: steps
1408 :
1409 : INTEGER :: i
1410 : REAL(KIND=dp) :: h, l, m
1411 :
1412 130 : l = lumo
1413 130 : h = homo
1414 :
1415 926 : DO i = 1, 200
1416 926 : IF (ABS(l) < threshold .AND. ABS(1 - h) < threshold) EXIT
1417 796 : m = 0.5_dp*(h + l)
1418 926 : IF (m > 0.5_dp) THEN
1419 412 : h = h**2
1420 412 : l = l**2
1421 : ELSE
1422 384 : h = 2*h - h**2
1423 384 : l = 2*l - l**2
1424 : END IF
1425 : END DO
1426 130 : steps = i - 1
1427 130 : END FUNCTION estimate_steps
1428 :
1429 : END MODULE dm_ls_scf_methods
|
