Line data Source code
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 12662 : 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 12662 : CALL timeset(routineN, handle)
83 :
84 : ! get a useful output_unit
85 12662 : logger => cp_get_default_logger()
86 12662 : IF (logger%para_env%is_source()) THEN
87 6331 : 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 12662 : IF (ls_scf_env%has_unit_metric) THEN
94 14 : CALL dbcsr_set(ls_scf_env%matrix_s, 0.0_dp)
95 14 : CALL dbcsr_add_on_diag(ls_scf_env%matrix_s, 1.0_dp)
96 : ELSE
97 12648 : CALL matrix_qs_to_ls(ls_scf_env%matrix_s, matrix_s, ls_scf_env%ls_mstruct, covariant=.TRUE.)
98 : END IF
99 :
100 12662 : CALL dbcsr_filter(ls_scf_env%matrix_s, ls_scf_env%eps_filter)
101 :
102 : ! needs a preconditioner for S
103 12662 : 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 290 : matrix_type=dbcsr_type_no_symmetry)
106 : CALL dbcsr_create(ls_scf_env%matrix_bs_sqrt_inv, template=ls_scf_env%matrix_s, &
107 290 : 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 290 : ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos)
113 : END IF
114 :
115 : ! precondition S
116 12662 : IF (ls_scf_env%has_s_preconditioner) THEN
117 : CALL apply_matrix_preconditioner(ls_scf_env%matrix_s, "forward", &
118 290 : 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 12662 : 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 12644 : matrix_type=dbcsr_type_no_symmetry)
126 : CALL dbcsr_create(ls_scf_env%matrix_s_sqrt_inv, template=ls_scf_env%matrix_s, &
127 12644 : matrix_type=dbcsr_type_no_symmetry)
128 :
129 12652 : 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 8 : 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 12636 : ls_scf_env%eps_lanczos, ls_scf_env%max_iter_lanczos)
141 : CASE DEFAULT
142 12644 : CPABORT("Unknown sqrt method.")
143 : END SELECT
144 :
145 12644 : IF (ls_scf_env%check_s_inv) THEN
146 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
147 0 : matrix_type=dbcsr_type_no_symmetry)
148 : CALL dbcsr_create(matrix_tmp2, template=ls_scf_env%matrix_s, &
149 0 : 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.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.0_dp, matrix_tmp2, filter_eps=ls_scf_env%eps_filter)
156 :
157 0 : frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp2)
158 0 : CALL dbcsr_add_on_diag(matrix_tmp2, -1.0_dp)
159 0 : frob_matrix = dbcsr_frobenius_norm(matrix_tmp2)
160 0 : IF (unit_nr > 0) THEN
161 0 : WRITE (unit_nr, *) "Error for (inv(sqrt(S))*S*inv(sqrt(S))-I)", frob_matrix/frob_matrix_base
162 : END IF
163 :
164 0 : CALL dbcsr_release(matrix_tmp1)
165 0 : CALL dbcsr_release(matrix_tmp2)
166 : END IF
167 : END IF
168 :
169 : ! compute the inverse of S
170 12662 : 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 12618 : matrix_type=dbcsr_type_no_symmetry)
173 12618 : IF (.NOT. ls_scf_env%use_s_sqrt) THEN
174 2 : 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 12616 : 0.0_dp, ls_scf_env%matrix_s_inv, filter_eps=ls_scf_env%eps_filter)
178 : END IF
179 12618 : IF (ls_scf_env%check_s_inv) THEN
180 : CALL dbcsr_create(matrix_tmp1, template=ls_scf_env%matrix_s, &
181 0 : 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.0_dp, matrix_tmp1, filter_eps=ls_scf_env%eps_filter)
184 0 : frob_matrix_base = dbcsr_frobenius_norm(matrix_tmp1)
185 0 : CALL dbcsr_add_on_diag(matrix_tmp1, -1.0_dp)
186 0 : frob_matrix = dbcsr_frobenius_norm(matrix_tmp1)
187 0 : IF (unit_nr > 0) THEN
188 0 : WRITE (unit_nr, *) "Error for (inv(S)*S-I)", frob_matrix/frob_matrix_base
189 : END IF
190 0 : CALL dbcsr_release(matrix_tmp1)
191 : END IF
192 : END IF
193 :
194 12662 : CALL timestop(handle)
195 12662 : 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 290 : 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 290 : REAL(dp), DIMENSION(:, :), POINTER :: block_dp
230 : TYPE(dbcsr_iterator_type) :: iter
231 : TYPE(dbcsr_type) :: matrix_bs
232 :
233 290 : 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 290 : CALL dbcsr_create(matrix_bs, template=matrix_s)
239 :
240 580 : SELECT CASE (preconditioner_type)
241 : CASE (ls_s_preconditioner_none)
242 : CASE (ls_s_preconditioner_atomic, ls_s_preconditioner_molecular)
243 290 : CALL dbcsr_iterator_start(iter, matrix_s)
244 27467 : DO WHILE (dbcsr_iterator_blocks_left(iter))
245 27177 : 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 27177 : block_needed = .FALSE.
253 :
254 27177 : IF (iblock_row == iblock_col) THEN
255 : block_needed = .TRUE.
256 : ELSE
257 24403 : IF (preconditioner_type == ls_s_preconditioner_molecular .AND. &
258 : ls_mstruct%cluster_type == ls_cluster_atomic) THEN
259 4756 : 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 290 : IF (block_needed) THEN
265 3962 : CALL dbcsr_put_block(matrix=matrix_bs, row=iblock_row, col=iblock_col, block=block_dp)
266 : END IF
267 :
268 : END DO
269 580 : CALL dbcsr_iterator_stop(iter)
270 : END SELECT
271 :
272 290 : CALL dbcsr_finalize(matrix_bs)
273 :
274 290 : SELECT CASE (preconditioner_type)
275 : CASE (ls_s_preconditioner_none)
276 : ! for now make it a simple identity matrix
277 0 : CALL dbcsr_copy(matrix_bs_sqrt, matrix_bs)
278 0 : CALL dbcsr_set(matrix_bs_sqrt, 0.0_dp)
279 0 : CALL dbcsr_add_on_diag(matrix_bs_sqrt, 1.0_dp)
280 :
281 : ! for now make it a simple identity matrix
282 0 : CALL dbcsr_copy(matrix_bs_sqrt_inv, matrix_bs)
283 0 : CALL dbcsr_set(matrix_bs_sqrt_inv, 0.0_dp)
284 0 : CALL dbcsr_add_on_diag(matrix_bs_sqrt_inv, 1.0_dp)
285 : CASE (ls_s_preconditioner_atomic, ls_s_preconditioner_molecular)
286 290 : CALL dbcsr_copy(matrix_bs_sqrt, matrix_bs)
287 290 : 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 580 : eps_lanczos=eps_lanczos, max_iter_lanczos=max_iter_lanczos)
295 : END SELECT
296 :
297 290 : CALL dbcsr_release(matrix_bs)
298 :
299 290 : CALL timestop(handle)
300 :
301 290 : 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 3278 : 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 3278 : CALL timeset(routineN, handle)
327 3278 : CALL dbcsr_create(matrix_tmp, template=matrix, matrix_type=dbcsr_type_no_symmetry)
328 :
329 2880 : SELECT CASE (direction)
330 : CASE ("forward")
331 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_bs_sqrt_inv, &
332 2880 : 0.0_dp, matrix_tmp)
333 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_bs_sqrt_inv, matrix_tmp, &
334 2880 : 0.0_dp, matrix)
335 : CASE ("backward")
336 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_bs_sqrt, &
337 398 : 0.0_dp, matrix_tmp)
338 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_bs_sqrt, matrix_tmp, &
339 398 : 0.0_dp, matrix)
340 : CASE DEFAULT
341 3278 : CPABORT("")
342 : END SELECT
343 :
344 3278 : CALL dbcsr_release(matrix_tmp)
345 :
346 3278 : CALL timestop(handle)
347 :
348 3278 : 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 910 : 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 910 : CALL timeset(routineN, handle)
397 :
398 910 : logger => cp_get_default_logger()
399 910 : IF (logger%para_env%is_source()) THEN
400 455 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
401 : ELSE
402 : unit_nr = -1
403 : END IF
404 :
405 910 : do_sign_symmetric = .FALSE.
406 910 : IF (PRESENT(sign_symmetric)) do_sign_symmetric = sign_symmetric
407 :
408 910 : used_submatrix_sign_method = ls_scf_submatrix_sign_ns
409 910 : 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 910 : used_submatrix_sign_method .EQ. ls_scf_submatrix_sign_direct_muadj_lowmem))
414 :
415 4 : 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 4 : nelectron, matrix_s_sqrt_inv)
420 : ELSE
421 906 : increment = initial_increment
422 :
423 906 : has_mu_low = .FALSE.
424 906 : has_mu_high = .FALSE.
425 :
426 : ! bisect if both bounds are known, otherwise find the bounds with a linear search
427 1050 : DO iter = 1, 30
428 1050 : IF (has_mu_low .AND. has_mu_high) THEN
429 16 : mu = (mu_low + mu_high)/2
430 16 : 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 1050 : matrix_s_sqrt_inv)
437 1050 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,I2,1X,F13.9,1X,F15.9)') &
438 525 : "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 1050 : IF (ABS(trace - nelectron) < 0.5_dp .OR. fixed_mu) EXIT
443 :
444 2100 : IF (trace < nelectron) THEN
445 32 : mu_low = mu
446 32 : mu = mu + increment
447 32 : has_mu_low = .TRUE.
448 32 : increment = increment*2
449 : ELSE
450 112 : mu_high = mu
451 112 : mu = mu - increment
452 112 : has_mu_high = .TRUE.
453 112 : increment = increment*2
454 : END IF
455 : END DO
456 :
457 : END IF
458 :
459 910 : CALL timestop(handle)
460 :
461 910 : 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 1072 : 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 1072 : CALL timeset(routineN, handle)
504 :
505 1072 : logger => cp_get_default_logger()
506 1072 : IF (logger%para_env%is_source()) THEN
507 536 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
508 : ELSE
509 : unit_nr = -1
510 : END IF
511 :
512 1072 : CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
513 :
514 1072 : IF (sign_symmetric) THEN
515 :
516 4 : IF (.NOT. PRESENT(matrix_s_sqrt_inv)) &
517 0 : CPABORT("Argument matrix_s_sqrt_inv required if sign_symmetric is set")
518 :
519 4 : CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
520 4 : 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 4 : 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 4 : 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
525 4 : CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)
526 :
527 6 : SELECT CASE (sign_method)
528 : CASE (ls_scf_sign_ns)
529 2 : CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
530 : CASE (ls_scf_sign_proot)
531 0 : CALL matrix_sign_proot(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order)
532 : CASE (ls_scf_sign_submatrix)
533 2 : CALL matrix_sign_submatrix(matrix_sign, matrix_ssqrtinv_ks_ssqrtinv, threshold, sign_order, submatrix_sign_method)
534 : CASE DEFAULT
535 4 : CPABORT("Unkown sign method.")
536 : END SELECT
537 4 : CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
538 4 : CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv2)
539 :
540 : ELSE ! .NOT. sign_symmetric
541 : ! get inv(S)*H-I*mu
542 1068 : 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 1068 : 0.0_dp, matrix_sinv_ks, filter_eps=threshold)
545 1068 : CALL dbcsr_add_on_diag(matrix_sinv_ks, -mu)
546 :
547 : ! compute sign(inv(S)*H-I*mu)
548 2124 : SELECT CASE (sign_method)
549 : CASE (ls_scf_sign_ns)
550 1056 : CALL matrix_sign_Newton_Schulz(matrix_sign, matrix_sinv_ks, threshold, sign_order)
551 : CASE (ls_scf_sign_proot)
552 8 : CALL matrix_sign_proot(matrix_sign, matrix_sinv_ks, threshold, sign_order)
553 : CASE (ls_scf_sign_submatrix)
554 4 : CALL matrix_sign_submatrix(matrix_sign, matrix_sinv_ks, threshold, sign_order, submatrix_sign_method)
555 : CASE DEFAULT
556 1068 : CPABORT("Unkown sign method.")
557 : END SELECT
558 1068 : 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 1072 : CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
563 1072 : CALL dbcsr_copy(matrix_p_ud, matrix_sign)
564 1072 : CALL dbcsr_scale(matrix_p_ud, -0.5_dp)
565 1072 : CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp)
566 1072 : CALL dbcsr_release(matrix_sign)
567 :
568 : ! we now have PS, lets get its trace
569 1072 : CALL dbcsr_trace(matrix_p_ud, trace)
570 :
571 : ! we can also check it is idempotent PS*PS=PS
572 1072 : 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 1072 : 0.0_dp, matrix_tmp, filter_eps=threshold)
575 1072 : CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp)
576 1072 : frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
577 1072 : IF (unit_nr > 0) WRITE (unit_nr, '(T2,A,F20.12)') "Deviation from idempotency: ", frob_matrix
578 :
579 1072 : IF (sign_symmetric) THEN
580 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_p_ud, &
581 4 : 0.0_dp, matrix_tmp, filter_eps=threshold)
582 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_s_sqrt_inv, &
583 4 : 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 1068 : 0.0_dp, matrix_p, filter_eps=threshold)
589 : END IF
590 1072 : CALL dbcsr_release(matrix_p_ud)
591 1072 : CALL dbcsr_release(matrix_tmp)
592 :
593 1072 : CALL timestop(handle)
594 :
595 1072 : 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 4 : 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 4 : CALL timeset(routineN, handle)
638 :
639 4 : logger => cp_get_default_logger()
640 4 : IF (logger%para_env%is_source()) THEN
641 2 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
642 : ELSE
643 : unit_nr = -1
644 : END IF
645 :
646 4 : CALL dbcsr_create(matrix_sign, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
647 :
648 4 : CALL dbcsr_create(matrix_ssqrtinv_ks_ssqrtinv, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
649 4 : 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 4 : 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 4 : 0.0_dp, matrix_ssqrtinv_ks_ssqrtinv, filter_eps=threshold)
654 4 : CALL dbcsr_add_on_diag(matrix_ssqrtinv_ks_ssqrtinv, -mu)
655 :
656 8 : SELECT CASE (sign_method)
657 : CASE (ls_scf_sign_submatrix)
658 8 : 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 4 : submatrix_sign_method)
662 : CASE DEFAULT
663 4 : CPABORT("density_matrix_sign_internal_mu called with invalid submatrix sign method")
664 : END SELECT
665 : CASE DEFAULT
666 4 : CPABORT("density_matrix_sign_internal_mu called with invalid sign method.")
667 : END SELECT
668 4 : CALL dbcsr_release(matrix_ssqrtinv_ks_ssqrtinv)
669 4 : 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 4 : CALL dbcsr_create(matrix_p_ud, template=matrix_s, matrix_type=dbcsr_type_no_symmetry)
673 4 : CALL dbcsr_copy(matrix_p_ud, matrix_sign)
674 4 : CALL dbcsr_scale(matrix_p_ud, -0.5_dp)
675 4 : CALL dbcsr_add_on_diag(matrix_p_ud, 0.5_dp)
676 4 : CALL dbcsr_release(matrix_sign)
677 :
678 : ! we now have PS, lets get its trace
679 4 : CALL dbcsr_trace(matrix_p_ud, trace)
680 :
681 : ! we can also check it is idempotent PS*PS=PS
682 4 : 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 4 : 0.0_dp, matrix_tmp, filter_eps=threshold)
685 4 : CALL dbcsr_add(matrix_tmp, matrix_p_ud, 1.0_dp, -1.0_dp)
686 4 : frob_matrix = dbcsr_frobenius_norm(matrix_tmp)
687 4 : 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 4 : 0.0_dp, matrix_tmp, filter_eps=threshold)
691 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_s_sqrt_inv, &
692 4 : 0.0_dp, matrix_p, filter_eps=threshold)
693 4 : CALL dbcsr_release(matrix_p_ud)
694 4 : CALL dbcsr_release(matrix_tmp)
695 :
696 4 : CALL timestop(handle)
697 :
698 4 : 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 13710 : 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 13710 : 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 13710 : IF (nelectron == 0) THEN
753 0 : CALL dbcsr_set(matrix_p, 0.0_dp)
754 : RETURN
755 : END IF
756 :
757 13710 : CALL timeset(routineN, handle)
758 :
759 13710 : logger => cp_get_default_logger()
760 13710 : IF (logger%para_env%is_source()) THEN
761 6855 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
762 : ELSE
763 6855 : unit_nr = -1
764 : END IF
765 :
766 13710 : do_dyn_threshold = .FALSE.
767 13710 : IF (PRESENT(dynamic_threshold)) do_dyn_threshold = dynamic_threshold
768 :
769 13710 : 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 13710 : 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 13710 : 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 13710 : 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 13710 : 0.0_dp, matrix_x, filter_eps=threshold)
781 13710 : CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym)
782 :
783 13710 : CALL dbcsr_create(matrix_k0, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
784 13710 : CALL dbcsr_copy(matrix_k0, matrix_x_nosym)
785 :
786 : ! compute the deviation in the mixed matrix, as seen in the ortho basis
787 13710 : IF (do_dyn_threshold) THEN
788 24 : CPASSERT(PRESENT(matrix_ks_deviation))
789 24 : 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 24 : converged=arnoldi_converged)
792 24 : maxdev = MAX(ABS(maxev), ABS(minev))
793 24 : IF (unit_nr > 0) THEN
794 12 : WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", arnoldi_converged
795 12 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "change in mixed matrix: ", maxdev
796 12 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo + maxdev
797 12 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo - maxdev
798 12 : 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 24 : 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 13710 : converged=arnoldi_converged)
808 20565 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", &
809 13710 : min_eig, max_eig, " converged: ", arnoldi_converged
810 13710 : eps_max = max_eig
811 13710 : eps_min = min_eig
812 :
813 : ! scale KS matrix
814 13710 : IF (eps_max == eps_min) THEN
815 20 : CALL dbcsr_scale(matrix_x, 1.0_dp/eps_max)
816 : ELSE
817 13690 : CALL dbcsr_add_on_diag(matrix_x, -eps_max)
818 13690 : CALL dbcsr_scale(matrix_x, -1.0_dp/(eps_max - eps_min))
819 : END IF
820 :
821 13710 : current_threshold = threshold
822 13710 : IF (do_dyn_threshold) THEN
823 : ! scale bounds for HOMO/LUMO
824 24 : scaled_homo_bound = (eps_max - (e_homo + maxdev))/(eps_max - eps_min)
825 24 : scaled_lumo_bound = (eps_max - (e_lumo - maxdev))/(eps_max - eps_min)
826 : END IF
827 :
828 13710 : CALL dbcsr_create(matrix_xsq, template=matrix_ks, matrix_type="S")
829 :
830 13710 : CALL dbcsr_create(matrix_xidsq, template=matrix_ks, matrix_type="S")
831 :
832 13710 : CALL dbcsr_create(tmp_gx, template=matrix_ks, matrix_type="S")
833 :
834 13710 : ALLOCATE (gamma_values(max_iter))
835 :
836 70448 : DO i = 1, max_iter
837 70448 : t1 = m_walltime()
838 70448 : 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 70448 : filter_eps=current_threshold, flop=flop1)
844 :
845 : ! intermediate use matrix_xidsq to compute = X*X-X
846 70448 : CALL dbcsr_copy(matrix_xidsq, matrix_x)
847 70448 : CALL dbcsr_add(matrix_xidsq, matrix_xsq, -1.0_dp, 1.0_dp)
848 70448 : frob_id = dbcsr_frobenius_norm(matrix_xidsq)
849 70448 : 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 70448 : CALL dbcsr_copy(matrix_xidsq, matrix_x)
854 70448 : CALL dbcsr_add(matrix_xidsq, matrix_xsq, -2.0_dp, 1.0_dp)
855 70448 : CALL dbcsr_add_on_diag(matrix_xidsq, 1.0_dp)
856 :
857 : ! tmp_gx = 4X-3X*X
858 70448 : CALL dbcsr_copy(tmp_gx, matrix_x)
859 70448 : 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 70448 : CALL dbcsr_dot(matrix_xsq, matrix_xidsq, trace_gx)
864 70448 : 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 70448 : delta_n = nelectron - trace_fx
870 : ! condition: ABS(frob_id/frob_x) < SQRT(threshold) ...
871 70448 : IF (((frob_id*frob_id) < (threshold*frob_x*frob_x)) .AND. (ABS(delta_n) < 0.5_dp)) THEN
872 13710 : gam = 3.0_dp
873 56738 : ELSE IF (ABS(delta_n) < 1e-14_dp) THEN
874 0 : 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 56738 : gam = delta_n/MAX(trace_gx, ABS(delta_n)/100)
878 : END IF
879 70448 : 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 70448 : IF (do_dyn_threshold) THEN
887 : ! quantities used for dynamic thresholding, when the estimated gap is larger than zero
888 154 : xi = (scaled_homo_bound - scaled_lumo_bound)
889 154 : IF (xi > 0.0_dp) THEN
890 130 : mmin = 0.5*(scaled_homo_bound + scaled_lumo_bound)
891 130 : max_threshold = ABS(1 - 2*mmin)*xi
892 :
893 130 : scaled_homo_bound = evaluate_trs4_polynomial(scaled_homo_bound, gamma_values(i:), 1)
894 130 : scaled_lumo_bound = evaluate_trs4_polynomial(scaled_lumo_bound, gamma_values(i:), 1)
895 130 : estimated_steps = estimate_steps(scaled_homo_bound, scaled_lumo_bound, threshold)
896 :
897 130 : est_threshold = (threshold/(estimated_steps + i + 1))*xi/(1 + threshold/(estimated_steps + i + 1))
898 130 : est_threshold = MIN(max_threshold, est_threshold)
899 130 : IF (i > 1) est_threshold = MAX(est_threshold, 0.1_dp*current_threshold)
900 130 : current_threshold = est_threshold
901 : ELSE
902 24 : current_threshold = threshold
903 : END IF
904 : END IF
905 :
906 70448 : IF (gam > gamma_max) THEN
907 : ! Xn+1 = 2X-X*X
908 908 : CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp)
909 908 : CALL dbcsr_filter(matrix_x, current_threshold)
910 908 : branch = 1
911 69540 : ELSE IF (gam < gamma_min) THEN
912 : ! Xn+1 = X*X
913 2996 : CALL dbcsr_copy(matrix_x, matrix_xsq)
914 2996 : branch = 2
915 : ELSE
916 : ! Xn+1 = F(X) + gam*G(X)
917 66544 : 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 66544 : flop=flop2, filter_eps=current_threshold)
921 66544 : branch = 3
922 : END IF
923 :
924 70448 : occ_matrix = dbcsr_get_occupation(matrix_x)
925 70448 : t2 = m_walltime()
926 70448 : IF (unit_nr > 0) THEN
927 : WRITE (unit_nr, &
928 35224 : '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TRS4 it ", &
929 35224 : i, occ_matrix, ABS(trace_gx), t2 - t1, &
930 70448 : (flop1 + flop2)/(1.0E6_dp*MAX(t2 - t1, 0.001_dp)), current_threshold
931 35224 : CALL m_flush(unit_nr)
932 : END IF
933 :
934 70448 : IF (abnormal_value(trace_gx)) &
935 0 : 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 70448 : IF ((frob_id*frob_id) < (threshold*frob_x*frob_x) .AND. branch == 3 .AND. (ABS(delta_n) < 0.5_dp)) THEN
941 13710 : IF (PRESENT(converged)) converged = .TRUE.
942 : EXIT
943 : END IF
944 :
945 : END DO
946 :
947 13710 : occ_matrix = dbcsr_get_occupation(matrix_x)
948 13710 : 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 13710 : CALL dbcsr_release(tmp_gx)
952 13710 : CALL dbcsr_release(matrix_xsq)
953 13710 : 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 13710 : 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 13710 : 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 13710 : mu_a = 0.0_dp; mu_b = 1.0_dp;
964 13710 : mu_fa = evaluate_trs4_polynomial(mu_a, gamma_values, i - 1) - 0.5_dp
965 161818 : DO j = 1, 40
966 161818 : mu_c = 0.5*(mu_a + mu_b)
967 161818 : 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 161818 : 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 161818 : IF (mu_fc*mu_fa > 0) THEN
971 76252 : mu_a = mu_c
972 76252 : mu_fa = mu_fc
973 : ELSE
974 : mu_b = mu_c
975 : END IF
976 : END DO
977 13710 : mu = (eps_min - eps_max)*mu_c + eps_max
978 13710 : DEALLOCATE (gamma_values)
979 13710 : IF (unit_nr > 0) THEN
980 6855 : WRITE (unit_nr, '(T6,A,1X,F12.5)') 'Chemical potential (mu): ', mu
981 : END IF
982 13710 : e_mu = mu
983 :
984 13710 : IF (do_dyn_threshold) THEN
985 24 : CALL dbcsr_desymmetrize(matrix_x, matrix_x_nosym)
986 : CALL compute_homo_lumo(matrix_k0, matrix_x_nosym, eps_min, eps_max, &
987 24 : threshold, max_iter_lanczos, eps_lanczos, homo, lumo, unit_nr)
988 24 : e_homo = homo
989 24 : e_lumo = lumo
990 : END IF
991 :
992 13710 : CALL dbcsr_release(matrix_x)
993 13710 : CALL dbcsr_release(matrix_x_nosym)
994 13710 : CALL dbcsr_release(matrix_k0)
995 13710 : CALL timestop(handle)
996 :
997 27420 : 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 50 : 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 50 : 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 50 : CALL timeset(routineN, handle)
1042 :
1043 50 : logger => cp_get_default_logger()
1044 50 : IF (logger%para_env%is_source()) THEN
1045 25 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
1046 : ELSE
1047 25 : unit_nr = -1
1048 : END IF
1049 :
1050 50 : do_non_monotonic = .FALSE.
1051 50 : 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 50 : CALL dbcsr_create(matrix_x, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1055 50 : 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 50 : 0.0_dp, matrix_xsq, filter_eps=threshold)
1059 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_xsq, matrix_s_sqrt_inv, &
1060 50 : 0.0_dp, matrix_x, filter_eps=threshold)
1061 :
1062 50 : IF (unit_nr > 0) THEN
1063 25 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "HOMO upper bound: ", e_homo
1064 25 : WRITE (unit_nr, '(T6,A,1X,F12.5)') "LUMO lower bound: ", e_lumo
1065 25 : 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 50 : converged=converged)
1071 75 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,1X,2F12.5,1X,A,1X,L1)') "Est. extremal eigenvalues", &
1072 50 : min_eig, max_eig, " converged: ", converged
1073 :
1074 50 : eps_max = max_eig
1075 50 : eps_min = min_eig
1076 :
1077 : ! scale KS matrix
1078 50 : CALL dbcsr_scale(matrix_x, -1.0_dp)
1079 50 : CALL dbcsr_add_on_diag(matrix_x, eps_max)
1080 50 : CALL dbcsr_scale(matrix_x, 1/(eps_max - eps_min))
1081 :
1082 50 : CALL dbcsr_copy(matrix_xsq, matrix_x)
1083 :
1084 50 : CALL dbcsr_create(matrix_tmp, template=matrix_ks, matrix_type=dbcsr_type_no_symmetry)
1085 :
1086 50 : ALLOCATE (poly(max_iter))
1087 50 : ALLOCATE (nu(max_iter))
1088 50 : ALLOCATE (wu(max_iter))
1089 50 : ALLOCATE (alpha(max_iter))
1090 50 : ALLOCATE (X(4))
1091 50 : ALLOCATE (Y(4))
1092 50 : ALLOCATE (lambda(4))
1093 :
1094 : ! Controls over the non_monotonic bounds, First if low gap, bias slightly
1095 50 : beta = (eps_max - ABS(e_lumo))/(eps_max - eps_min)
1096 50 : betaB = (eps_max + ABS(e_homo))/(eps_max - eps_min)
1097 :
1098 50 : IF ((beta - betaB) < 0.005_dp) THEN
1099 50 : beta = beta - 0.002_dp
1100 50 : betaB = betaB + 0.002_dp
1101 : END IF
1102 : ! Check if input specifies to use monotonic bounds.
1103 50 : IF (.NOT. do_non_monotonic) THEN
1104 24 : beta = 0.0_dp
1105 24 : betaB = 1.0_dp
1106 : END IF
1107 : ! initial SCF cycle has no reliable estimate of homo/lumo, force monotinic bounds.
1108 50 : IF (e_homo == 0.0_dp) THEN
1109 10 : beta = 0.0_dp
1110 10 : BetaB = 1.0_dp
1111 : END IF
1112 :
1113 : ! init to take true branch first
1114 50 : trace_fx = nelectron
1115 50 : trace_gx = 0
1116 :
1117 802 : DO i = 1, max_iter
1118 802 : t1 = m_walltime()
1119 802 : flop1 = 0; flop2 = 0
1120 :
1121 802 : IF (ABS(trace_fx - nelectron) <= ABS(trace_gx - nelectron)) THEN
1122 : ! Xn+1 = (aX+ (1-a)I)^2
1123 406 : poly(i) = 1.0_dp
1124 406 : alpha(i) = 2.0_dp/(2.0_dp - beta)
1125 :
1126 406 : CALL dbcsr_scale(matrix_x, alpha(i))
1127 406 : 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 406 : filter_eps=threshold, flop=flop1)
1131 :
1132 : !save X for control variables
1133 406 : CALL dbcsr_copy(matrix_tmp, matrix_x)
1134 :
1135 406 : CALL dbcsr_copy(matrix_x, matrix_xsq)
1136 :
1137 406 : beta = (1.0_dp - alpha(i)) + alpha(i)*beta
1138 406 : beta = beta*beta
1139 406 : betaB = (1.0_dp - alpha(i)) + alpha(i)*betaB
1140 406 : betaB = betaB*betaB
1141 : ELSE
1142 : ! Xn+1 = 2aX-a^2*X^2
1143 396 : poly(i) = 0.0_dp
1144 396 : alpha(i) = 2.0_dp/(1.0_dp + betaB)
1145 :
1146 396 : 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 396 : filter_eps=threshold, flop=flop1)
1150 :
1151 : !save X for control variables
1152 396 : CALL dbcsr_copy(matrix_tmp, matrix_x)
1153 : !
1154 396 : CALL dbcsr_add(matrix_x, matrix_xsq, 2.0_dp, -1.0_dp)
1155 :
1156 396 : beta = alpha(i)*beta
1157 396 : beta = 2.0_dp*beta - beta*beta
1158 396 : betaB = alpha(i)*betaB
1159 396 : betaB = 2.0_dp*betaB - betaB*betaB
1160 :
1161 : END IF
1162 802 : occ_matrix = dbcsr_get_occupation(matrix_x)
1163 802 : t2 = m_walltime()
1164 802 : IF (unit_nr > 0) THEN
1165 : WRITE (unit_nr, &
1166 401 : '(T6,A,I3,1X,F10.8,E12.3,F12.3,F13.3,E12.3)') "TC2 it ", &
1167 401 : i, occ_matrix, t2 - t1, &
1168 802 : (flop1 + flop2)/(1.0E6_dp*(t2 - t1)), threshold
1169 401 : CALL m_flush(unit_nr)
1170 : END IF
1171 :
1172 : ! calculate control terms
1173 802 : CALL dbcsr_trace(matrix_xsq, trace_fx)
1174 :
1175 : ! intermediate use matrix_xsq compute X- X*X , temorarily use trace_gx
1176 802 : CALL dbcsr_add(matrix_xsq, matrix_tmp, -1.0_dp, 1.0_dp)
1177 802 : CALL dbcsr_trace(matrix_xsq, trace_gx)
1178 802 : nu(i) = dbcsr_frobenius_norm(matrix_xsq)
1179 802 : wu(i) = trace_gx
1180 :
1181 : ! intermediate use matrix_xsq to compute = 2X - X*X
1182 802 : CALL dbcsr_add(matrix_xsq, matrix_tmp, 1.0_dp, 1.0_dp)
1183 802 : CALL dbcsr_trace(matrix_xsq, trace_gx)
1184 : ! TC2 has quadratic convergence, using the frobeniums norm as an idempotency deviation test.
1185 2406 : IF (ABS(nu(i)) < (threshold)) EXIT
1186 : END DO
1187 :
1188 50 : occ_matrix = dbcsr_get_occupation(matrix_x)
1189 50 : 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 50 : 0.0_dp, matrix_tmp, filter_eps=threshold)
1194 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_s_sqrt_inv, matrix_tmp, &
1195 50 : 0.0_dp, matrix_p, filter_eps=threshold)
1196 :
1197 50 : CALL dbcsr_release(matrix_xsq)
1198 50 : CALL dbcsr_release(matrix_tmp)
1199 :
1200 : ! ALGO 3 from. SIAM DOI. 10.1137/130911585
1201 50 : X(1) = 1.0_dp
1202 50 : X(2) = 1.0_dp
1203 50 : X(3) = 0.0_dp
1204 50 : X(4) = 0.0_dp
1205 : gama = 6.0_dp - 4.0_dp*(SQRT(2.0_dp))
1206 50 : gama = gama - gama*gama
1207 504 : DO WHILE (nu(i) < gama)
1208 : ! safeguard against negative root, is skipping correct?
1209 454 : IF (wu(i) < 1.0e-14_dp) THEN
1210 108 : i = i - 1
1211 108 : CYCLE
1212 : END IF
1213 346 : IF ((1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)) < 0.0_dp) THEN
1214 18 : i = i - 1
1215 18 : CYCLE
1216 : END IF
1217 328 : Y(1) = 0.5_dp*(1.0_dp - SQRT(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)))
1218 328 : Y(2) = 0.5_dp*(1.0_dp - SQRT(1.0_dp - 4.0_dp*nu(i)))
1219 328 : Y(3) = 0.5_dp*(1.0_dp + SQRT(1.0_dp - 4.0_dp*nu(i)))
1220 328 : Y(4) = 0.5_dp*(1.0_dp + SQRT(1.0_dp - 4.0_dp*nu(i)*nu(i)/wu(i)))
1221 1640 : Y(:) = MIN(1.0_dp, MAX(0.0_dp, Y(:)))
1222 3978 : DO j = i, 1, -1
1223 3978 : IF (poly(j) == 1.0_dp) THEN
1224 9110 : DO k = 1, 4
1225 7288 : Y(k) = SQRT(Y(k))
1226 9110 : Y(k) = (Y(k) - 1.0_dp + alpha(j))/alpha(j)
1227 : END DO ! end K
1228 : ELSE
1229 9140 : DO k = 1, 4
1230 7312 : Y(k) = 1.0_dp - SQRT(1.0_dp - Y(k))
1231 9140 : Y(k) = Y(k)/alpha(j)
1232 : END DO ! end K
1233 : END IF ! end poly
1234 : END DO ! end j
1235 328 : X(1) = MIN(X(1), Y(1))
1236 328 : X(2) = MIN(X(2), Y(2))
1237 328 : X(3) = MAX(X(3), Y(3))
1238 328 : X(4) = MAX(X(4), Y(4))
1239 328 : i = i - 1
1240 : END DO ! end i
1241 : ! lambda 1,2,3,4 are:: out lumo, in lumo, in homo, out homo
1242 250 : DO k = 1, 4
1243 250 : 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 50 : e_homo = lambda(4)
1247 50 : e_lumo = lambda(1)
1248 50 : IF (unit_nr > 0) WRITE (unit_nr, '(T6,A,3E12.4)') "outer homo/lumo/gap", e_homo, e_lumo, (e_lumo - e_homo)
1249 50 : 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 50 : DEALLOCATE (poly)
1252 50 : DEALLOCATE (nu)
1253 50 : DEALLOCATE (wu)
1254 50 : DEALLOCATE (alpha)
1255 50 : DEALLOCATE (X)
1256 50 : DEALLOCATE (Y)
1257 50 : DEALLOCATE (lambda)
1258 :
1259 50 : CALL dbcsr_release(matrix_x)
1260 50 : CALL timestop(handle)
1261 :
1262 100 : 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 132 : 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 44 : CALL dbcsr_create(tmp1, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1296 :
1297 44 : CALL dbcsr_create(tmp2, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1298 :
1299 44 : CALL dbcsr_create(tmp3, template=matrix_k, matrix_type=dbcsr_type_no_symmetry)
1300 :
1301 44 : shift1 = -eps_min
1302 44 : 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 44 : CALL dbcsr_copy(tmp2, matrix_k)
1306 44 : CALL dbcsr_add_on_diag(tmp2, shift1)
1307 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p, tmp2, &
1308 44 : 0.0_dp, tmp1, filter_eps=threshold)
1309 : CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, &
1310 44 : threshold=eps_lanczos, max_iter=max_iter_lanczos)
1311 44 : homo = max_eig - shift1
1312 44 : IF (unit_nr > 0) THEN
1313 22 : 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 44 : CALL dbcsr_copy(tmp3, matrix_p)
1318 44 : CALL dbcsr_scale(tmp3, -1.0_dp)
1319 44 : CALL dbcsr_add_on_diag(tmp3, 1.0_dp) !tmp3 = 1-P
1320 44 : CALL dbcsr_copy(tmp2, matrix_k)
1321 44 : CALL dbcsr_add_on_diag(tmp2, -shift2)
1322 : CALL dbcsr_multiply("N", "N", -1.0_dp, tmp3, tmp2, &
1323 44 : 0.0_dp, tmp1, filter_eps=threshold)
1324 : CALL arnoldi_extremal(tmp1, max_eig, min_eig, converged=converged, &
1325 44 : threshold=eps_lanczos, max_iter=max_iter_lanczos)
1326 44 : lumo = -max_eig + shift2
1327 :
1328 44 : IF (unit_nr > 0) THEN
1329 22 : WRITE (unit_nr, '(T6,A,1X,L12)') "Lanczos converged: ", converged
1330 22 : WRITE (unit_nr, '(T6,A,1X,3F12.5)') 'HOMO/LUMO/gap', homo, lumo, lumo - homo
1331 : END IF
1332 44 : CALL dbcsr_release(tmp1)
1333 44 : CALL dbcsr_release(tmp2)
1334 44 : CALL dbcsr_release(tmp3)
1335 :
1336 44 : END SUBROUTINE compute_homo_lumo
1337 :
1338 : ! **************************************************************************************************
1339 : !> \brief ...
1340 : !> \param x ...
1341 : !> \param gamma_values ...
1342 : !> \param i ...
1343 : !> \return ...
1344 : ! **************************************************************************************************
1345 175788 : 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 175788 : xr = x
1356 1362012 : DO k = 1, i
1357 1362012 : IF (gamma_values(k) > gam_max) THEN
1358 18990 : xr = 2*xr - xr**2
1359 1167234 : ELSE IF (gamma_values(k) < gam_min) THEN
1360 62690 : xr = xr**2
1361 : ELSE
1362 1104544 : xr = (xr*xr)*(4*xr - 3*xr*xr) + gamma_values(k)*xr*xr*((1 - xr)**2)
1363 : END IF
1364 : END DO
1365 175788 : END FUNCTION evaluate_trs4_polynomial
1366 :
1367 : ! **************************************************************************************************
1368 : !> \brief ...
1369 : !> \param homo ...
1370 : !> \param lumo ...
1371 : !> \param threshold ...
1372 : !> \return ...
1373 : ! **************************************************************************************************
1374 130 : 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 130 : l = lumo
1382 130 : h = homo
1383 :
1384 926 : DO i = 1, 200
1385 926 : IF (ABS(l) < threshold .AND. ABS(1 - h) < threshold) EXIT
1386 796 : m = 0.5_dp*(h + l)
1387 926 : IF (m > 0.5_dp) THEN
1388 412 : h = h**2
1389 412 : l = l**2
1390 : ELSE
1391 384 : h = 2*h - h**2
1392 384 : l = 2*l - l**2
1393 : END IF
1394 : END DO
1395 130 : steps = i - 1
1396 130 : END FUNCTION estimate_steps
1397 :
1398 : END MODULE dm_ls_scf_methods
|
