Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Cayley transformation methods
10 : !> \par History
11 : !> 2011.06 created [Rustam Z Khaliullin]
12 : !> \author Rustam Z Khaliullin
13 : ! **************************************************************************************************
14 : MODULE ct_methods
15 : USE cp_dbcsr_api, ONLY: &
16 : dbcsr_add, dbcsr_add_on_diag, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, &
17 : dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_dot, dbcsr_filter, dbcsr_finalize, &
18 : dbcsr_get_diag, dbcsr_get_info, dbcsr_get_stored_coordinates, dbcsr_iterator_blocks_left, &
19 : dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
20 : dbcsr_multiply, dbcsr_nblkcols_total, dbcsr_nblkrows_total, dbcsr_release, &
21 : dbcsr_reserve_block2d, dbcsr_scale, dbcsr_set, dbcsr_set_diag, dbcsr_transposed, &
22 : dbcsr_type, dbcsr_type_no_symmetry, dbcsr_work_create
23 : USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,&
24 : cp_dbcsr_cholesky_invert
25 : USE cp_dbcsr_contrib, ONLY: dbcsr_frobenius_norm,&
26 : dbcsr_hadamard_product,&
27 : dbcsr_maxabs
28 : USE cp_dbcsr_diag, ONLY: cp_dbcsr_syevd
29 : USE cp_log_handling, ONLY: cp_get_default_logger,&
30 : cp_logger_get_default_unit_nr,&
31 : cp_logger_type
32 : USE ct_types, ONLY: ct_step_env_type
33 : USE input_constants, ONLY: &
34 : cg_dai_yuan, cg_fletcher, cg_fletcher_reeves, cg_hager_zhang, cg_hestenes_stiefel, &
35 : cg_liu_storey, cg_polak_ribiere, cg_zero, tensor_orthogonal, tensor_up_down
36 : USE iterate_matrix, ONLY: matrix_sqrt_Newton_Schulz
37 : USE kinds, ONLY: dp
38 : USE machine, ONLY: m_walltime
39 : #include "./base/base_uses.f90"
40 :
41 : IMPLICIT NONE
42 :
43 : PRIVATE
44 :
45 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ct_methods'
46 :
47 : ! Public subroutines
48 : PUBLIC :: ct_step_execute, analytic_line_search, diagonalize_diagonal_blocks
49 :
50 : CONTAINS
51 :
52 : ! **************************************************************************************************
53 : !> \brief Performs Cayley transformation
54 : !> \param cts_env ...
55 : !> \par History
56 : !> 2011.06 created [Rustam Z Khaliullin]
57 : !> \author Rustam Z Khaliullin
58 : ! **************************************************************************************************
59 0 : SUBROUTINE ct_step_execute(cts_env)
60 :
61 : TYPE(ct_step_env_type) :: cts_env
62 :
63 : CHARACTER(len=*), PARAMETER :: routineN = 'ct_step_execute'
64 :
65 : INTEGER :: handle, n, preconditioner_type, unit_nr
66 : REAL(KIND=dp) :: gap_estimate, safety_margin
67 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: evals
68 : TYPE(cp_logger_type), POINTER :: logger
69 : TYPE(dbcsr_type) :: matrix_pp, matrix_pq, matrix_qp, &
70 : matrix_qp_save, matrix_qq, oo1, &
71 : oo1_sqrt, oo1_sqrt_inv, t_corr, tmp1, &
72 : u_pp, u_qq
73 :
74 : !TYPE(dbcsr_type) :: rst_x1, rst_x2
75 : !REAL(KIND=dp) :: ener_tmp
76 : !TYPE(dbcsr_iterator_type) :: iter
77 : !INTEGER :: iblock_row,iblock_col,&
78 : ! iblock_row_size,iblock_col_size
79 : !REAL(KIND=dp), DIMENSION(:,:), POINTER :: data_p
80 :
81 0 : CALL timeset(routineN, handle)
82 :
83 0 : logger => cp_get_default_logger()
84 0 : IF (logger%para_env%is_source()) THEN
85 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
86 : ELSE
87 : unit_nr = -1
88 : END IF
89 :
90 : ! check if all input is in place and flags are consistent
91 0 : IF (cts_env%update_q .AND. (.NOT. cts_env%update_p)) THEN
92 0 : CPABORT("q-update is possible only with p-update")
93 : END IF
94 :
95 0 : IF (cts_env%tensor_type .EQ. tensor_up_down) THEN
96 0 : CPABORT("riccati is not implemented for biorthogonal basis")
97 : END IF
98 :
99 0 : IF (.NOT. ASSOCIATED(cts_env%matrix_ks)) THEN
100 0 : CPABORT("KS matrix is not associated")
101 : END IF
102 :
103 0 : IF (cts_env%use_virt_orbs .AND. (.NOT. cts_env%use_occ_orbs)) THEN
104 0 : CPABORT("virtual orbs can be used only with occupied orbs")
105 : END IF
106 :
107 0 : IF (cts_env%use_occ_orbs) THEN
108 0 : IF (.NOT. ASSOCIATED(cts_env%matrix_t)) THEN
109 0 : CPABORT("T matrix is not associated")
110 : END IF
111 0 : IF (.NOT. ASSOCIATED(cts_env%matrix_qp_template)) THEN
112 0 : CPABORT("QP template is not associated")
113 : END IF
114 0 : IF (.NOT. ASSOCIATED(cts_env%matrix_pq_template)) THEN
115 0 : CPABORT("PQ template is not associated")
116 : END IF
117 : END IF
118 :
119 0 : IF (cts_env%use_virt_orbs) THEN
120 0 : IF (.NOT. ASSOCIATED(cts_env%matrix_v)) THEN
121 0 : CPABORT("V matrix is not associated")
122 : END IF
123 : ELSE
124 0 : IF (.NOT. ASSOCIATED(cts_env%matrix_p)) THEN
125 0 : CPABORT("P matrix is not associated")
126 : END IF
127 : END IF
128 :
129 0 : IF (cts_env%tensor_type .NE. tensor_up_down .AND. &
130 : cts_env%tensor_type .NE. tensor_orthogonal) THEN
131 0 : CPABORT("illegal tensor flag")
132 : END IF
133 :
134 : ! start real calculations
135 0 : IF (cts_env%use_occ_orbs) THEN
136 :
137 : ! create matrices for various ks blocks
138 : CALL dbcsr_create(matrix_pp, &
139 : template=cts_env%p_index_up, &
140 0 : matrix_type=dbcsr_type_no_symmetry)
141 : CALL dbcsr_create(matrix_qp, &
142 : template=cts_env%matrix_qp_template, &
143 0 : matrix_type=dbcsr_type_no_symmetry)
144 : CALL dbcsr_create(matrix_qq, &
145 : template=cts_env%q_index_up, &
146 0 : matrix_type=dbcsr_type_no_symmetry)
147 : CALL dbcsr_create(matrix_pq, &
148 : template=cts_env%matrix_pq_template, &
149 0 : matrix_type=dbcsr_type_no_symmetry)
150 :
151 : ! create the residue matrix
152 : CALL dbcsr_create(cts_env%matrix_res, &
153 0 : template=cts_env%matrix_qp_template)
154 :
155 : CALL assemble_ks_qp_blocks(cts_env%matrix_ks, &
156 : cts_env%matrix_p, &
157 : cts_env%matrix_t, &
158 : cts_env%matrix_v, &
159 : cts_env%q_index_down, &
160 : cts_env%p_index_up, &
161 : cts_env%q_index_up, &
162 : matrix_pp, &
163 : matrix_qq, &
164 : matrix_qp, &
165 : matrix_pq, &
166 : cts_env%tensor_type, &
167 : cts_env%use_virt_orbs, &
168 0 : cts_env%eps_filter)
169 :
170 : ! create a matrix of single-excitation amplitudes
171 : CALL dbcsr_create(cts_env%matrix_x, &
172 0 : template=cts_env%matrix_qp_template)
173 0 : IF (ASSOCIATED(cts_env%matrix_x_guess)) THEN
174 : CALL dbcsr_copy(cts_env%matrix_x, &
175 0 : cts_env%matrix_x_guess)
176 0 : IF (cts_env%tensor_type .EQ. tensor_orthogonal) THEN
177 : ! bring x from contravariant-covariant representation
178 : ! to the orthogonal/cholesky representation
179 : ! use res as temporary storage
180 : CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_down, &
181 : cts_env%matrix_x, 0.0_dp, cts_env%matrix_res, &
182 0 : filter_eps=cts_env%eps_filter)
183 : CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_res, &
184 : cts_env%p_index_up, 0.0_dp, &
185 : cts_env%matrix_x, &
186 0 : filter_eps=cts_env%eps_filter)
187 : END IF
188 : ELSE
189 : ! set amplitudes to zero
190 0 : CALL dbcsr_set(cts_env%matrix_x, 0.0_dp)
191 : END IF
192 :
193 : !SELECT CASE (cts_env%preconditioner_type)
194 : !CASE (prec_eigenvector_blocks,prec_eigenvector_full)
195 0 : preconditioner_type = 1
196 0 : safety_margin = 2.0_dp
197 0 : gap_estimate = 0.0001_dp
198 : SELECT CASE (preconditioner_type)
199 : CASE (1, 2)
200 : !RZK-warning diagonalization works only with orthogonal tensor!!!
201 : ! find a better basis by diagonalizing diagonal blocks
202 : ! first pp
203 : CALL dbcsr_create(u_pp, template=matrix_pp, &
204 0 : matrix_type=dbcsr_type_no_symmetry)
205 : !IF (cts_env%preconditioner_type.eq.prec_eigenvector_full) THEN
206 : IF (.TRUE.) THEN
207 0 : CALL dbcsr_get_info(matrix_pp, nfullrows_total=n)
208 0 : ALLOCATE (evals(n))
209 : CALL cp_dbcsr_syevd(matrix_pp, u_pp, evals, &
210 0 : cts_env%para_env, cts_env%blacs_env)
211 0 : DEALLOCATE (evals)
212 : ELSE
213 : CALL diagonalize_diagonal_blocks(matrix_pp, u_pp)
214 : END IF
215 : ! and now qq
216 : CALL dbcsr_create(u_qq, template=matrix_qq, &
217 0 : matrix_type=dbcsr_type_no_symmetry)
218 : !IF (cts_env%preconditioner_type.eq.prec_eigenvector_full) THEN
219 : IF (.TRUE.) THEN
220 0 : CALL dbcsr_get_info(matrix_qq, nfullrows_total=n)
221 0 : ALLOCATE (evals(n))
222 : CALL cp_dbcsr_syevd(matrix_qq, u_qq, evals, &
223 0 : cts_env%para_env, cts_env%blacs_env)
224 0 : DEALLOCATE (evals)
225 : ELSE
226 : CALL diagonalize_diagonal_blocks(matrix_qq, u_qq)
227 : END IF
228 :
229 : ! apply the transformation to all matrices
230 : CALL matrix_forward_transform(matrix_pp, u_pp, u_pp, &
231 0 : cts_env%eps_filter)
232 : CALL matrix_forward_transform(matrix_qq, u_qq, u_qq, &
233 0 : cts_env%eps_filter)
234 : CALL matrix_forward_transform(matrix_qp, u_qq, u_pp, &
235 0 : cts_env%eps_filter)
236 : CALL matrix_forward_transform(matrix_pq, u_pp, u_qq, &
237 0 : cts_env%eps_filter)
238 : CALL matrix_forward_transform(cts_env%matrix_x, u_qq, u_pp, &
239 0 : cts_env%eps_filter)
240 :
241 0 : IF (cts_env%max_iter .GE. 0) THEN
242 :
243 : CALL solve_riccati_equation( &
244 : pp=matrix_pp, &
245 : qq=matrix_qq, &
246 : qp=matrix_qp, &
247 : pq=matrix_pq, &
248 : x=cts_env%matrix_x, &
249 : res=cts_env%matrix_res, &
250 : neglect_quadratic_term=cts_env%neglect_quadratic_term, &
251 : conjugator=cts_env%conjugator, &
252 : max_iter=cts_env%max_iter, &
253 : eps_convergence=cts_env%eps_convergence, &
254 : eps_filter=cts_env%eps_filter, &
255 0 : converged=cts_env%converged)
256 :
257 0 : IF (cts_env%converged) THEN
258 : !IF (unit_nr>0) THEN
259 : ! WRITE(unit_nr,*)
260 : ! WRITE(unit_nr,'(T6,A)') &
261 : ! "RICCATI equations solved"
262 : ! CALL m_flush(unit_nr)
263 : !ENDIF
264 : ELSE
265 0 : CPABORT("RICCATI: CG algorithm has NOT converged")
266 : END IF
267 :
268 : END IF
269 :
270 0 : IF (cts_env%calculate_energy_corr) THEN
271 :
272 0 : CALL dbcsr_dot(matrix_qp, cts_env%matrix_x, cts_env%energy_correction)
273 :
274 : END IF
275 :
276 0 : CALL dbcsr_release(matrix_pp)
277 0 : CALL dbcsr_release(matrix_qp)
278 0 : CALL dbcsr_release(matrix_qq)
279 0 : CALL dbcsr_release(matrix_pq)
280 :
281 : ! back-transform to the original basis
282 : CALL matrix_backward_transform(cts_env%matrix_x, u_qq, &
283 0 : u_pp, cts_env%eps_filter)
284 :
285 0 : CALL dbcsr_release(u_qq)
286 0 : CALL dbcsr_release(u_pp)
287 :
288 : !CASE (prec_cholesky_inverse)
289 : CASE (3)
290 :
291 : ! RZK-warning implemented only for orthogonal tensors!!!
292 : ! generalization to up_down should be easy
293 : CALL dbcsr_create(u_pp, template=matrix_pp, &
294 : matrix_type=dbcsr_type_no_symmetry)
295 : CALL dbcsr_copy(u_pp, matrix_pp)
296 : CALL dbcsr_scale(u_pp, -1.0_dp)
297 : CALL dbcsr_add_on_diag(u_pp, &
298 : ABS(safety_margin*gap_estimate))
299 : CALL cp_dbcsr_cholesky_decompose(u_pp, &
300 : para_env=cts_env%para_env, &
301 : blacs_env=cts_env%blacs_env)
302 : CALL cp_dbcsr_cholesky_invert(u_pp, &
303 : para_env=cts_env%para_env, &
304 : blacs_env=cts_env%blacs_env, &
305 : uplo_to_full=.TRUE.)
306 : !CALL dbcsr_scale(u_pp,-1.0_dp)
307 :
308 : CALL dbcsr_create(u_qq, template=matrix_qq, &
309 : matrix_type=dbcsr_type_no_symmetry)
310 : CALL dbcsr_copy(u_qq, matrix_qq)
311 : CALL dbcsr_add_on_diag(u_qq, &
312 : ABS(safety_margin*gap_estimate))
313 : CALL cp_dbcsr_cholesky_decompose(u_qq, &
314 : para_env=cts_env%para_env, &
315 : blacs_env=cts_env%blacs_env)
316 : CALL cp_dbcsr_cholesky_invert(u_qq, &
317 : para_env=cts_env%para_env, &
318 : blacs_env=cts_env%blacs_env, &
319 : uplo_to_full=.TRUE.)
320 :
321 : ! transform all riccati matrices (left-right preconditioner)
322 : CALL dbcsr_create(tmp1, template=matrix_qq, &
323 : matrix_type=dbcsr_type_no_symmetry)
324 : CALL dbcsr_multiply("N", "N", 1.0_dp, u_qq, &
325 : matrix_qq, 0.0_dp, tmp1, &
326 : filter_eps=cts_env%eps_filter)
327 : CALL dbcsr_copy(matrix_qq, tmp1)
328 : CALL dbcsr_release(tmp1)
329 :
330 : CALL dbcsr_create(tmp1, template=matrix_pp, &
331 : matrix_type=dbcsr_type_no_symmetry)
332 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_pp, &
333 : u_pp, 0.0_dp, tmp1, &
334 : filter_eps=cts_env%eps_filter)
335 : CALL dbcsr_copy(matrix_pp, tmp1)
336 : CALL dbcsr_release(tmp1)
337 :
338 : CALL dbcsr_create(matrix_qp_save, template=matrix_qp, &
339 : matrix_type=dbcsr_type_no_symmetry)
340 : CALL dbcsr_copy(matrix_qp_save, matrix_qp)
341 :
342 : CALL dbcsr_create(tmp1, template=matrix_qp, &
343 : matrix_type=dbcsr_type_no_symmetry)
344 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
345 : u_pp, 0.0_dp, tmp1, &
346 : filter_eps=cts_env%eps_filter)
347 : CALL dbcsr_multiply("N", "N", 1.0_dp, u_qq, tmp1, &
348 : 0.0_dp, matrix_qp, &
349 : filter_eps=cts_env%eps_filter)
350 : CALL dbcsr_release(tmp1)
351 : !CALL dbcsr_print(matrix_qq)
352 : !CALL dbcsr_print(matrix_qp)
353 : !CALL dbcsr_print(matrix_pp)
354 :
355 : IF (cts_env%max_iter .GE. 0) THEN
356 :
357 : CALL solve_riccati_equation( &
358 : pp=matrix_pp, &
359 : qq=matrix_qq, &
360 : qp=matrix_qp, &
361 : pq=matrix_pq, &
362 : oo=u_pp, &
363 : vv=u_qq, &
364 : x=cts_env%matrix_x, &
365 : res=cts_env%matrix_res, &
366 : neglect_quadratic_term=cts_env%neglect_quadratic_term, &
367 : conjugator=cts_env%conjugator, &
368 : max_iter=cts_env%max_iter, &
369 : eps_convergence=cts_env%eps_convergence, &
370 : eps_filter=cts_env%eps_filter, &
371 : converged=cts_env%converged)
372 :
373 : IF (cts_env%converged) THEN
374 : !IF (unit_nr>0) THEN
375 : ! WRITE(unit_nr,*)
376 : ! WRITE(unit_nr,'(T6,A)') &
377 : ! "RICCATI equations solved"
378 : ! CALL m_flush(unit_nr)
379 : !ENDIF
380 : ELSE
381 : CPABORT("RICCATI: CG algorithm has NOT converged")
382 : END IF
383 :
384 : END IF
385 :
386 : IF (cts_env%calculate_energy_corr) THEN
387 :
388 : CALL dbcsr_dot(matrix_qp_save, cts_env%matrix_x, cts_env%energy_correction)
389 :
390 : END IF
391 : CALL dbcsr_release(matrix_qp_save)
392 :
393 : CALL dbcsr_release(matrix_pp)
394 : CALL dbcsr_release(matrix_qp)
395 : CALL dbcsr_release(matrix_qq)
396 : CALL dbcsr_release(matrix_pq)
397 :
398 : CALL dbcsr_release(u_qq)
399 : CALL dbcsr_release(u_pp)
400 :
401 : CASE DEFAULT
402 : CPABORT("illegal preconditioner type")
403 : END SELECT ! preconditioner type
404 :
405 0 : IF (cts_env%update_p) THEN
406 :
407 0 : IF (cts_env%tensor_type .EQ. tensor_up_down) THEN
408 0 : CPABORT("orbital update is NYI for this tensor type")
409 : END IF
410 :
411 : ! transform occupied orbitals
412 : ! in a way that preserves the overlap metric
413 : CALL dbcsr_create(oo1, &
414 : template=cts_env%p_index_up, &
415 0 : matrix_type=dbcsr_type_no_symmetry)
416 : CALL dbcsr_create(oo1_sqrt_inv, &
417 0 : template=oo1)
418 : CALL dbcsr_create(oo1_sqrt, &
419 0 : template=oo1)
420 :
421 : ! Compute (1+tr(X).X)^(-1/2)_up_down
422 : CALL dbcsr_multiply("T", "N", 1.0_dp, cts_env%matrix_x, &
423 : cts_env%matrix_x, 0.0_dp, oo1, &
424 0 : filter_eps=cts_env%eps_filter)
425 0 : CALL dbcsr_add_on_diag(oo1, 1.0_dp)
426 : CALL matrix_sqrt_Newton_Schulz(oo1_sqrt, &
427 : oo1_sqrt_inv, &
428 : oo1, &
429 : !if cholesky is used then sqrt
430 : !guess cannot be provided
431 : !matrix_sqrt_inv_guess=cts_env%p_index_up,&
432 : !matrix_sqrt_guess=cts_env%p_index_down,&
433 : threshold=cts_env%eps_filter, &
434 : order=cts_env%order_lanczos, &
435 : eps_lanczos=cts_env%eps_lancsoz, &
436 0 : max_iter_lanczos=cts_env%max_iter_lanczos)
437 : CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%p_index_up, &
438 : oo1_sqrt_inv, 0.0_dp, oo1, &
439 0 : filter_eps=cts_env%eps_filter)
440 : CALL dbcsr_multiply("N", "N", 1.0_dp, oo1, &
441 : cts_env%p_index_down, 0.0_dp, oo1_sqrt, &
442 0 : filter_eps=cts_env%eps_filter)
443 0 : CALL dbcsr_release(oo1)
444 0 : CALL dbcsr_release(oo1_sqrt_inv)
445 :
446 : ! bring x to contravariant-covariant representation now
447 : CALL dbcsr_create(matrix_qp, &
448 : template=cts_env%matrix_qp_template, &
449 0 : matrix_type=dbcsr_type_no_symmetry)
450 : CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_up, &
451 : cts_env%matrix_x, 0.0_dp, matrix_qp, &
452 0 : filter_eps=cts_env%eps_filter)
453 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
454 : cts_env%p_index_down, 0.0_dp, &
455 : cts_env%matrix_x, &
456 0 : filter_eps=cts_env%eps_filter)
457 0 : CALL dbcsr_release(matrix_qp)
458 :
459 : ! update T=T+X or T=T+V.X (whichever is appropriate)
460 0 : CALL dbcsr_create(t_corr, template=cts_env%matrix_t)
461 0 : IF (cts_env%use_virt_orbs) THEN
462 : CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_v, &
463 : cts_env%matrix_x, 0.0_dp, t_corr, &
464 0 : filter_eps=cts_env%eps_filter)
465 : CALL dbcsr_add(cts_env%matrix_t, t_corr, &
466 0 : 1.0_dp, 1.0_dp)
467 : ELSE
468 : CALL dbcsr_add(cts_env%matrix_t, cts_env%matrix_x, &
469 0 : 1.0_dp, 1.0_dp)
470 : END IF
471 : ! adjust T so the metric is preserved: T=(T+X).(1+tr(X).X)^(-1/2)
472 : CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%matrix_t, oo1_sqrt, &
473 0 : 0.0_dp, t_corr, filter_eps=cts_env%eps_filter)
474 0 : CALL dbcsr_copy(cts_env%matrix_t, t_corr)
475 :
476 0 : CALL dbcsr_release(t_corr)
477 0 : CALL dbcsr_release(oo1_sqrt)
478 :
479 : ELSE ! do not update p
480 :
481 0 : IF (cts_env%tensor_type .EQ. tensor_orthogonal) THEN
482 : ! bring x to contravariant-covariant representation
483 : CALL dbcsr_create(matrix_qp, &
484 : template=cts_env%matrix_qp_template, &
485 0 : matrix_type=dbcsr_type_no_symmetry)
486 : CALL dbcsr_multiply("N", "N", 1.0_dp, cts_env%q_index_up, &
487 : cts_env%matrix_x, 0.0_dp, matrix_qp, &
488 0 : filter_eps=cts_env%eps_filter)
489 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_qp, &
490 : cts_env%p_index_down, 0.0_dp, &
491 : cts_env%matrix_x, &
492 0 : filter_eps=cts_env%eps_filter)
493 0 : CALL dbcsr_release(matrix_qp)
494 : END IF
495 :
496 : END IF
497 :
498 : ELSE
499 0 : CPABORT("illegal occ option")
500 : END IF
501 :
502 0 : CALL timestop(handle)
503 :
504 0 : END SUBROUTINE ct_step_execute
505 :
506 : ! **************************************************************************************************
507 : !> \brief computes oo, ov, vo, and vv blocks of the ks matrix
508 : !> \param ks ...
509 : !> \param p ...
510 : !> \param t ...
511 : !> \param v ...
512 : !> \param q_index_down ...
513 : !> \param p_index_up ...
514 : !> \param q_index_up ...
515 : !> \param pp ...
516 : !> \param qq ...
517 : !> \param qp ...
518 : !> \param pq ...
519 : !> \param tensor_type ...
520 : !> \param use_virt_orbs ...
521 : !> \param eps_filter ...
522 : !> \par History
523 : !> 2011.06 created [Rustam Z Khaliullin]
524 : !> \author Rustam Z Khaliullin
525 : ! **************************************************************************************************
526 0 : SUBROUTINE assemble_ks_qp_blocks(ks, p, t, v, q_index_down, &
527 : p_index_up, q_index_up, pp, qq, qp, pq, tensor_type, use_virt_orbs, eps_filter)
528 :
529 : TYPE(dbcsr_type), INTENT(IN) :: ks, p, t, v, q_index_down, p_index_up, &
530 : q_index_up
531 : TYPE(dbcsr_type), INTENT(OUT) :: pp, qq, qp, pq
532 : INTEGER, INTENT(IN) :: tensor_type
533 : LOGICAL, INTENT(IN) :: use_virt_orbs
534 : REAL(KIND=dp), INTENT(IN) :: eps_filter
535 :
536 : CHARACTER(len=*), PARAMETER :: routineN = 'assemble_ks_qp_blocks'
537 :
538 : INTEGER :: handle
539 : LOGICAL :: library_fixed
540 : TYPE(dbcsr_type) :: kst, ksv, no, on, oo, q_index_up_nosym, &
541 : sp, spf, t_or, v_or
542 :
543 0 : CALL timeset(routineN, handle)
544 :
545 0 : IF (use_virt_orbs) THEN
546 :
547 : ! orthogonalize the orbitals
548 0 : CALL dbcsr_create(t_or, template=t)
549 0 : CALL dbcsr_create(v_or, template=v)
550 : CALL dbcsr_multiply("N", "N", 1.0_dp, t, p_index_up, &
551 0 : 0.0_dp, t_or, filter_eps=eps_filter)
552 : CALL dbcsr_multiply("N", "N", 1.0_dp, v, q_index_up, &
553 0 : 0.0_dp, v_or, filter_eps=eps_filter)
554 :
555 : ! KS.T
556 0 : CALL dbcsr_create(kst, template=t)
557 : CALL dbcsr_multiply("N", "N", 1.0_dp, ks, t_or, &
558 0 : 0.0_dp, kst, filter_eps=eps_filter)
559 : ! pp=tr(T)*KS.T
560 : CALL dbcsr_multiply("T", "N", 1.0_dp, t_or, kst, &
561 0 : 0.0_dp, pp, filter_eps=eps_filter)
562 : ! qp=tr(V)*KS.T
563 : CALL dbcsr_multiply("T", "N", 1.0_dp, v_or, kst, &
564 0 : 0.0_dp, qp, filter_eps=eps_filter)
565 0 : CALL dbcsr_release(kst)
566 :
567 : ! KS.V
568 0 : CALL dbcsr_create(ksv, template=v)
569 : CALL dbcsr_multiply("N", "N", 1.0_dp, ks, v_or, &
570 0 : 0.0_dp, ksv, filter_eps=eps_filter)
571 : ! tr(T)*KS.V
572 : CALL dbcsr_multiply("T", "N", 1.0_dp, t_or, ksv, &
573 0 : 0.0_dp, pq, filter_eps=eps_filter)
574 : ! tr(V)*KS.V
575 : CALL dbcsr_multiply("T", "N", 1.0_dp, v_or, ksv, &
576 0 : 0.0_dp, qq, filter_eps=eps_filter)
577 0 : CALL dbcsr_release(ksv)
578 :
579 0 : CALL dbcsr_release(t_or)
580 0 : CALL dbcsr_release(v_or)
581 :
582 : ELSE ! no virtuals, use projected AOs
583 :
584 : ! THIS PROCEDURE HAS NOT BEEN UPDATED FOR CHOLESKY p/q_index_up/down
585 : CALL dbcsr_create(sp, template=q_index_down, &
586 0 : matrix_type=dbcsr_type_no_symmetry)
587 : CALL dbcsr_create(spf, template=q_index_down, &
588 0 : matrix_type=dbcsr_type_no_symmetry)
589 :
590 : ! qp=KS*T
591 : CALL dbcsr_multiply("N", "N", 1.0_dp, ks, t, 0.0_dp, qp, &
592 0 : filter_eps=eps_filter)
593 : ! pp=tr(T)*KS.T
594 : CALL dbcsr_multiply("T", "N", 1.0_dp, t, qp, 0.0_dp, pp, &
595 0 : filter_eps=eps_filter)
596 : ! sp=-S_*P
597 : CALL dbcsr_multiply("N", "N", -1.0_dp, q_index_down, p, 0.0_dp, sp, &
598 0 : filter_eps=eps_filter)
599 :
600 : ! sp=1/S^-S_.P
601 0 : SELECT CASE (tensor_type)
602 : CASE (tensor_up_down)
603 0 : CALL dbcsr_add_on_diag(sp, 1.0_dp)
604 : CASE (tensor_orthogonal)
605 : CALL dbcsr_create(q_index_up_nosym, template=q_index_up, &
606 0 : matrix_type=dbcsr_type_no_symmetry)
607 0 : CALL dbcsr_desymmetrize(q_index_up, q_index_up_nosym)
608 0 : CALL dbcsr_add(sp, q_index_up_nosym, 1.0_dp, 1.0_dp)
609 0 : CALL dbcsr_release(q_index_up_nosym)
610 : END SELECT
611 :
612 : ! spf=(1/S^-S_.P)*KS
613 : CALL dbcsr_multiply("N", "N", 1.0_dp, sp, ks, 0.0_dp, spf, &
614 0 : filter_eps=eps_filter)
615 :
616 : ! qp=spf*T
617 : CALL dbcsr_multiply("N", "N", 1.0_dp, spf, t, 0.0_dp, qp, &
618 0 : filter_eps=eps_filter)
619 :
620 0 : SELECT CASE (tensor_type)
621 : CASE (tensor_up_down)
622 : ! pq=tr(qp)
623 0 : CALL dbcsr_transposed(pq, qp, transpose_distribution=.FALSE.)
624 : CASE (tensor_orthogonal)
625 : ! pq=sig^.tr(qp)
626 : CALL dbcsr_multiply("N", "T", 1.0_dp, p_index_up, qp, 0.0_dp, pq, &
627 0 : filter_eps=eps_filter)
628 0 : library_fixed = .FALSE.
629 0 : IF (library_fixed) THEN
630 : CALL dbcsr_transposed(qp, pq, transpose_distribution=.FALSE.)
631 : ELSE
632 : CALL dbcsr_create(no, template=qp, &
633 0 : matrix_type=dbcsr_type_no_symmetry)
634 : CALL dbcsr_multiply("N", "N", 1.0_dp, qp, p_index_up, 0.0_dp, no, &
635 0 : filter_eps=eps_filter)
636 0 : CALL dbcsr_copy(qp, no)
637 0 : CALL dbcsr_release(no)
638 : END IF
639 : END SELECT
640 :
641 : ! qq=spf*tr(sp)
642 : CALL dbcsr_multiply("N", "T", 1.0_dp, spf, sp, 0.0_dp, qq, &
643 0 : filter_eps=eps_filter)
644 :
645 0 : SELECT CASE (tensor_type)
646 : CASE (tensor_up_down)
647 :
648 : CALL dbcsr_create(oo, template=pp, &
649 0 : matrix_type=dbcsr_type_no_symmetry)
650 : CALL dbcsr_create(no, template=qp, &
651 0 : matrix_type=dbcsr_type_no_symmetry)
652 :
653 : ! first index up
654 : CALL dbcsr_multiply("N", "N", 1.0_dp, q_index_up, qq, 0.0_dp, spf, &
655 0 : filter_eps=eps_filter)
656 0 : CALL dbcsr_copy(qq, spf)
657 : CALL dbcsr_multiply("N", "N", 1.0_dp, q_index_up, qp, 0.0_dp, no, &
658 0 : filter_eps=eps_filter)
659 0 : CALL dbcsr_copy(qp, no)
660 : CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pp, 0.0_dp, oo, &
661 0 : filter_eps=eps_filter)
662 0 : CALL dbcsr_copy(pp, oo)
663 : CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pq, 0.0_dp, on, &
664 0 : filter_eps=eps_filter)
665 0 : CALL dbcsr_copy(pq, on)
666 :
667 0 : CALL dbcsr_release(no)
668 0 : CALL dbcsr_release(oo)
669 :
670 : CASE (tensor_orthogonal)
671 :
672 : CALL dbcsr_create(oo, template=pp, &
673 0 : matrix_type=dbcsr_type_no_symmetry)
674 :
675 : ! both indeces up in the pp block
676 : CALL dbcsr_multiply("N", "N", 1.0_dp, p_index_up, pp, 0.0_dp, oo, &
677 0 : filter_eps=eps_filter)
678 : CALL dbcsr_multiply("N", "N", 1.0_dp, oo, p_index_up, 0.0_dp, pp, &
679 0 : filter_eps=eps_filter)
680 :
681 0 : CALL dbcsr_release(oo)
682 :
683 : END SELECT
684 :
685 0 : CALL dbcsr_release(sp)
686 0 : CALL dbcsr_release(spf)
687 :
688 : END IF
689 :
690 0 : CALL timestop(handle)
691 :
692 0 : END SUBROUTINE assemble_ks_qp_blocks
693 :
694 : ! **************************************************************************************************
695 : !> \brief Solves the generalized Riccati or Sylvester eqation
696 : !> using the preconditioned conjugate gradient algorithm
697 : !> qp + qq.x.oo - vv.x.pp - vv.x.pq.x.oo = 0 [oo and vv are optional]
698 : !> qp + qq.x - x.pp - x.pq.x = 0
699 : !> \param pp ...
700 : !> \param qq ...
701 : !> \param qp ...
702 : !> \param pq ...
703 : !> \param oo ...
704 : !> \param vv ...
705 : !> \param x ...
706 : !> \param res ...
707 : !> \param neglect_quadratic_term ...
708 : !> \param conjugator ...
709 : !> \param max_iter ...
710 : !> \param eps_convergence ...
711 : !> \param eps_filter ...
712 : !> \param converged ...
713 : !> \par History
714 : !> 2011.06 created [Rustam Z Khaliullin]
715 : !> 2011.11 generalized [Rustam Z Khaliullin]
716 : !> \author Rustam Z Khaliullin
717 : ! **************************************************************************************************
718 0 : RECURSIVE SUBROUTINE solve_riccati_equation(pp, qq, qp, pq, oo, vv, x, res, &
719 : neglect_quadratic_term, &
720 : conjugator, max_iter, eps_convergence, eps_filter, &
721 : converged)
722 :
723 : TYPE(dbcsr_type), INTENT(IN) :: pp, qq
724 : TYPE(dbcsr_type), INTENT(INOUT) :: qp
725 : TYPE(dbcsr_type), INTENT(IN) :: pq
726 : TYPE(dbcsr_type), INTENT(IN), OPTIONAL :: oo, vv
727 : TYPE(dbcsr_type), INTENT(INOUT) :: x
728 : TYPE(dbcsr_type), INTENT(OUT) :: res
729 : LOGICAL, INTENT(IN) :: neglect_quadratic_term
730 : INTEGER, INTENT(IN) :: conjugator, max_iter
731 : REAL(KIND=dp), INTENT(IN) :: eps_convergence, eps_filter
732 : LOGICAL, INTENT(OUT) :: converged
733 :
734 : CHARACTER(len=*), PARAMETER :: routineN = 'solve_riccati_equation'
735 :
736 : INTEGER :: handle, istep, iteration, nsteps, &
737 : unit_nr, update_prec_freq
738 : LOGICAL :: prepare_to_exit, present_oo, present_vv, &
739 : quadratic_term, restart_conjugator
740 : REAL(KIND=dp) :: best_norm, best_step_size, beta, c0, c1, &
741 : c2, c3, denom, kappa, numer, &
742 : obj_function, t1, t2, tau
743 : REAL(KIND=dp), DIMENSION(3) :: step_size
744 : TYPE(cp_logger_type), POINTER :: logger
745 : TYPE(dbcsr_type) :: aux1, aux2, grad, m, n, oo1, oo2, prec, &
746 : res_trial, step, step_oo, vv_step
747 :
748 : !TYPE(dbcsr_type) :: qqqq, pppp, zero_pq, zero_qp
749 :
750 0 : CALL timeset(routineN, handle)
751 :
752 0 : logger => cp_get_default_logger()
753 0 : IF (logger%para_env%is_source()) THEN
754 0 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
755 : ELSE
756 : unit_nr = -1
757 : END IF
758 :
759 0 : t1 = m_walltime()
760 :
761 : !IF (level.gt.5) THEN
762 : ! CPErrorMessage(cp_failure_level,routineP,"recursion level is too high")
763 : ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
764 : !ENDIF
765 : !IF (unit_nr>0) THEN
766 : ! WRITE(unit_nr,*) &
767 : ! "========== LEVEL ",level,"=========="
768 : !ENDIF
769 : !CALL dbcsr_print(qq)
770 : !CALL dbcsr_print(pp)
771 : !CALL dbcsr_print(qp)
772 : !!CALL dbcsr_print(pq)
773 : !IF (unit_nr>0) THEN
774 : ! WRITE(unit_nr,*) &
775 : ! "====== END LEVEL ",level,"=========="
776 : !ENDIF
777 :
778 0 : quadratic_term = .NOT. neglect_quadratic_term
779 0 : present_oo = PRESENT(oo)
780 0 : present_vv = PRESENT(vv)
781 :
782 : ! create aux1 matrix and init
783 0 : CALL dbcsr_create(aux1, template=pp)
784 0 : CALL dbcsr_copy(aux1, pp)
785 0 : CALL dbcsr_scale(aux1, -1.0_dp)
786 :
787 : ! create aux2 matrix and init
788 0 : CALL dbcsr_create(aux2, template=qq)
789 0 : CALL dbcsr_copy(aux2, qq)
790 :
791 : ! create the gradient matrix and init
792 0 : CALL dbcsr_create(grad, template=x)
793 0 : CALL dbcsr_set(grad, 0.0_dp)
794 :
795 : ! create a preconditioner
796 : ! RZK-warning how to apply it to up_down tensor?
797 0 : CALL dbcsr_create(prec, template=x)
798 : !CALL create_preconditioner(prec,aux1,aux2,qp,res,tensor_type,eps_filter)
799 : !CALL dbcsr_set(prec,1.0_dp)
800 :
801 : ! create the step matrix and init
802 0 : CALL dbcsr_create(step, template=x)
803 : !CALL dbcsr_hadamard_product(prec,grad,step)
804 : !CALL dbcsr_scale(step,-1.0_dp)
805 :
806 0 : CALL dbcsr_create(n, template=x)
807 0 : CALL dbcsr_create(m, template=x)
808 0 : CALL dbcsr_create(oo1, template=pp)
809 0 : CALL dbcsr_create(oo2, template=pp)
810 0 : CALL dbcsr_create(res_trial, template=res)
811 0 : CALL dbcsr_create(vv_step, template=res)
812 0 : CALL dbcsr_create(step_oo, template=res)
813 :
814 : ! start conjugate gradient iterations
815 0 : iteration = 0
816 0 : converged = .FALSE.
817 0 : prepare_to_exit = .FALSE.
818 0 : beta = 0.0_dp
819 0 : best_step_size = 0.0_dp
820 0 : best_norm = 1.0E+100_dp
821 : !ecorr=0.0_dp
822 : !change_ecorr=0.0_dp
823 0 : restart_conjugator = .FALSE.
824 0 : update_prec_freq = 20
825 : DO
826 :
827 : ! (re)-compute the residuals
828 0 : IF (iteration .EQ. 0) THEN
829 0 : CALL dbcsr_copy(res, qp)
830 0 : IF (present_oo) THEN
831 : CALL dbcsr_multiply("N", "N", +1.0_dp, qq, x, 0.0_dp, res_trial, &
832 0 : filter_eps=eps_filter)
833 : CALL dbcsr_multiply("N", "N", +1.0_dp, res_trial, oo, 1.0_dp, res, &
834 0 : filter_eps=eps_filter)
835 : ELSE
836 : CALL dbcsr_multiply("N", "N", +1.0_dp, qq, x, 1.0_dp, res, &
837 0 : filter_eps=eps_filter)
838 : END IF
839 0 : IF (present_vv) THEN
840 : CALL dbcsr_multiply("N", "N", -1.0_dp, x, pp, 0.0_dp, res_trial, &
841 0 : filter_eps=eps_filter)
842 : CALL dbcsr_multiply("N", "N", +1.0_dp, vv, res_trial, 1.0_dp, res, &
843 0 : filter_eps=eps_filter)
844 : ELSE
845 : CALL dbcsr_multiply("N", "N", -1.0_dp, x, pp, 1.0_dp, res, &
846 0 : filter_eps=eps_filter)
847 : END IF
848 0 : IF (quadratic_term) THEN
849 0 : IF (present_oo) THEN
850 : CALL dbcsr_multiply("N", "N", +1.0_dp, pq, x, 0.0_dp, oo1, &
851 0 : filter_eps=eps_filter)
852 : CALL dbcsr_multiply("N", "N", +1.0_dp, oo1, oo, 0.0_dp, oo2, &
853 0 : filter_eps=eps_filter)
854 : ELSE
855 : CALL dbcsr_multiply("N", "N", +1.0_dp, pq, x, 0.0_dp, oo2, &
856 0 : filter_eps=eps_filter)
857 : END IF
858 0 : IF (present_vv) THEN
859 : CALL dbcsr_multiply("N", "N", -1.0_dp, x, oo2, 0.0_dp, res_trial, &
860 0 : filter_eps=eps_filter)
861 : CALL dbcsr_multiply("N", "N", +1.0_dp, vv, res_trial, 1.0_dp, res, &
862 0 : filter_eps=eps_filter)
863 : ELSE
864 : CALL dbcsr_multiply("N", "N", -1.0_dp, x, oo2, 1.0_dp, res, &
865 0 : filter_eps=eps_filter)
866 : END IF
867 : END IF
868 0 : best_norm = dbcsr_maxabs(res)
869 : ELSE
870 0 : CALL dbcsr_add(res, m, 1.0_dp, best_step_size)
871 0 : CALL dbcsr_add(res, n, 1.0_dp, -best_step_size*best_step_size)
872 0 : CALL dbcsr_filter(res, eps_filter)
873 : END IF
874 :
875 : ! check convergence and other exit criteria
876 0 : converged = (best_norm .LT. eps_convergence)
877 0 : IF (converged .OR. (iteration .GE. max_iter)) THEN
878 : prepare_to_exit = .TRUE.
879 : END IF
880 :
881 0 : IF (.NOT. prepare_to_exit) THEN
882 :
883 : ! update aux1=-pp-pq.x.oo and aux2=qq-vv.x.pq
884 0 : IF (quadratic_term) THEN
885 0 : IF (iteration .EQ. 0) THEN
886 0 : IF (present_oo) THEN
887 : CALL dbcsr_multiply("N", "N", -1.0_dp, pq, x, 0.0_dp, oo1, &
888 0 : filter_eps=eps_filter)
889 : CALL dbcsr_multiply("N", "N", +1.0_dp, oo1, oo, 1.0_dp, aux1, &
890 0 : filter_eps=eps_filter)
891 : ELSE
892 : CALL dbcsr_multiply("N", "N", -1.0_dp, pq, x, 1.0_dp, aux1, &
893 0 : filter_eps=eps_filter)
894 : END IF
895 0 : IF (present_vv) THEN
896 : CALL dbcsr_multiply("N", "N", -1.0_dp, vv, x, 0.0_dp, res_trial, &
897 0 : filter_eps=eps_filter)
898 : CALL dbcsr_multiply("N", "N", +1.0_dp, res_trial, pq, 1.0_dp, aux2, &
899 0 : filter_eps=eps_filter)
900 : ELSE
901 : CALL dbcsr_multiply("N", "N", -1.0_dp, x, pq, 1.0_dp, aux2, &
902 0 : filter_eps=eps_filter)
903 : END IF
904 : ELSE
905 0 : IF (present_oo) THEN
906 : CALL dbcsr_multiply("N", "N", -best_step_size, pq, step_oo, 1.0_dp, aux1, &
907 0 : filter_eps=eps_filter)
908 : ELSE
909 : CALL dbcsr_multiply("N", "N", -best_step_size, pq, step, 1.0_dp, aux1, &
910 0 : filter_eps=eps_filter)
911 : END IF
912 0 : IF (present_vv) THEN
913 : CALL dbcsr_multiply("N", "N", -best_step_size, vv_step, pq, 1.0_dp, aux2, &
914 0 : filter_eps=eps_filter)
915 : ELSE
916 : CALL dbcsr_multiply("N", "N", -best_step_size, step, pq, 1.0_dp, aux2, &
917 0 : filter_eps=eps_filter)
918 : END IF
919 : END IF
920 : END IF
921 :
922 : ! recompute the gradient, do not update it yet
923 : ! use m matrix as a temporary storage
924 : ! grad=t(vv).res.t(aux1)+t(aux2).res.t(oo)
925 0 : IF (present_vv) THEN
926 : CALL dbcsr_multiply("N", "T", 1.0_dp, res, aux1, 0.0_dp, res_trial, &
927 0 : filter_eps=eps_filter)
928 : CALL dbcsr_multiply("T", "N", 1.0_dp, vv, res_trial, 0.0_dp, m, &
929 0 : filter_eps=eps_filter)
930 : ELSE
931 : CALL dbcsr_multiply("N", "T", 1.0_dp, res, aux1, 0.0_dp, m, &
932 0 : filter_eps=eps_filter)
933 : END IF
934 0 : IF (present_oo) THEN
935 : CALL dbcsr_multiply("T", "N", 1.0_dp, aux1, res, 0.0_dp, res_trial, &
936 0 : filter_eps=eps_filter)
937 : CALL dbcsr_multiply("N", "T", 1.0_dp, res_trial, oo, 1.0_dp, m, &
938 0 : filter_eps=eps_filter)
939 : ELSE
940 : CALL dbcsr_multiply("T", "N", 1.0_dp, aux2, res, 1.0_dp, m, &
941 0 : filter_eps=eps_filter)
942 : END IF
943 :
944 : ! compute preconditioner
945 : !IF (iteration.eq.0.OR.(mod(iteration,update_prec_freq).eq.0)) THEN
946 0 : IF (iteration .EQ. 0) THEN
947 0 : CALL create_preconditioner(prec, aux1, aux2, eps_filter)
948 : !restart_conjugator=.TRUE.
949 : !CALL dbcsr_set(prec,1.0_dp)
950 : !CALL dbcsr_print(prec)
951 : END IF
952 :
953 : ! compute the conjugation coefficient - beta
954 0 : IF ((iteration .EQ. 0) .OR. restart_conjugator) THEN
955 0 : beta = 0.0_dp
956 : ELSE
957 0 : restart_conjugator = .FALSE.
958 0 : SELECT CASE (conjugator)
959 : CASE (cg_hestenes_stiefel)
960 0 : CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
961 0 : CALL dbcsr_hadamard_product(prec, grad, n)
962 0 : CALL dbcsr_dot(n, m, numer)
963 0 : CALL dbcsr_dot(grad, step, denom)
964 0 : beta = numer/denom
965 : CASE (cg_fletcher_reeves)
966 0 : CALL dbcsr_hadamard_product(prec, grad, n)
967 0 : CALL dbcsr_dot(grad, n, denom)
968 0 : CALL dbcsr_hadamard_product(prec, m, n)
969 0 : CALL dbcsr_dot(m, n, numer)
970 0 : beta = numer/denom
971 : CASE (cg_polak_ribiere)
972 0 : CALL dbcsr_hadamard_product(prec, grad, n)
973 0 : CALL dbcsr_dot(grad, n, denom)
974 0 : CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
975 0 : CALL dbcsr_hadamard_product(prec, grad, n)
976 0 : CALL dbcsr_dot(n, m, numer)
977 0 : beta = numer/denom
978 : CASE (cg_fletcher)
979 0 : CALL dbcsr_hadamard_product(prec, m, n)
980 0 : CALL dbcsr_dot(m, n, numer)
981 0 : CALL dbcsr_dot(grad, step, denom)
982 0 : beta = -1.0_dp*numer/denom
983 : CASE (cg_liu_storey)
984 0 : CALL dbcsr_dot(grad, step, denom)
985 0 : CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
986 0 : CALL dbcsr_hadamard_product(prec, grad, n)
987 0 : CALL dbcsr_dot(n, m, numer)
988 0 : beta = -1.0_dp*numer/denom
989 : CASE (cg_dai_yuan)
990 0 : CALL dbcsr_hadamard_product(prec, m, n)
991 0 : CALL dbcsr_dot(m, n, numer)
992 0 : CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
993 0 : CALL dbcsr_dot(grad, step, denom)
994 0 : beta = numer/denom
995 : CASE (cg_hager_zhang)
996 0 : CALL dbcsr_add(grad, m, -1.0_dp, 1.0_dp)
997 0 : CALL dbcsr_dot(grad, step, denom)
998 0 : CALL dbcsr_hadamard_product(prec, grad, n)
999 0 : CALL dbcsr_dot(n, grad, numer)
1000 0 : kappa = 2.0_dp*numer/denom
1001 0 : CALL dbcsr_dot(n, m, numer)
1002 0 : tau = numer/denom
1003 0 : CALL dbcsr_dot(step, m, numer)
1004 0 : beta = tau - kappa*numer/denom
1005 : CASE (cg_zero)
1006 0 : beta = 0.0_dp
1007 : CASE DEFAULT
1008 0 : CPABORT("illegal conjugator")
1009 : END SELECT
1010 : END IF ! iteration.eq.0
1011 :
1012 : ! move the current gradient to its storage
1013 0 : CALL dbcsr_copy(grad, m)
1014 :
1015 : ! precondition new gradient (use m as tmp storage)
1016 0 : CALL dbcsr_hadamard_product(prec, grad, m)
1017 0 : CALL dbcsr_filter(m, eps_filter)
1018 :
1019 : ! recompute the step direction
1020 0 : CALL dbcsr_add(step, m, beta, -1.0_dp)
1021 0 : CALL dbcsr_filter(step, eps_filter)
1022 :
1023 : !! ALTERNATIVE METHOD TO OBTAIN THE STEP FROM THE GRADIENT
1024 : !CALL dbcsr_init(qqqq)
1025 : !CALL dbcsr_create(qqqq,template=qq)
1026 : !CALL dbcsr_init(pppp)
1027 : !CALL dbcsr_create(pppp,template=pp)
1028 : !CALL dbcsr_init(zero_pq)
1029 : !CALL dbcsr_create(zero_pq,template=pq)
1030 : !CALL dbcsr_init(zero_qp)
1031 : !CALL dbcsr_create(zero_qp,template=qp)
1032 : !CALL dbcsr_multiply("T","N",1.0_dp,aux2,aux2,0.0_dp,qqqq,&
1033 : ! filter_eps=eps_filter)
1034 : !CALL dbcsr_multiply("N","T",-1.0_dp,aux1,aux1,0.0_dp,pppp,&
1035 : ! filter_eps=eps_filter)
1036 : !CALL dbcsr_set(zero_qp,0.0_dp)
1037 : !CALL dbcsr_set(zero_pq,0.0_dp)
1038 : !CALL solve_riccati_equation(pppp,qqqq,grad,zero_pq,zero_qp,zero_qp,&
1039 : ! .TRUE.,tensor_type,&
1040 : ! conjugator,max_iter,eps_convergence,eps_filter,&
1041 : ! converged,level+1)
1042 : !CALL dbcsr_release(qqqq)
1043 : !CALL dbcsr_release(pppp)
1044 : !CALL dbcsr_release(zero_qp)
1045 : !CALL dbcsr_release(zero_pq)
1046 :
1047 : ! calculate the optimal step size
1048 : ! m=step.aux1+aux2.step
1049 0 : IF (present_vv) THEN
1050 : CALL dbcsr_multiply("N", "N", 1.0_dp, vv, step, 0.0_dp, vv_step, &
1051 0 : filter_eps=eps_filter)
1052 : CALL dbcsr_multiply("N", "N", 1.0_dp, vv_step, aux1, 0.0_dp, m, &
1053 0 : filter_eps=eps_filter)
1054 : ELSE
1055 : CALL dbcsr_multiply("N", "N", 1.0_dp, step, aux1, 0.0_dp, m, &
1056 0 : filter_eps=eps_filter)
1057 : END IF
1058 0 : IF (present_oo) THEN
1059 : CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo, 0.0_dp, step_oo, &
1060 0 : filter_eps=eps_filter)
1061 : CALL dbcsr_multiply("N", "N", 1.0_dp, aux2, step_oo, 1.0_dp, m, &
1062 0 : filter_eps=eps_filter)
1063 : ELSE
1064 : CALL dbcsr_multiply("N", "N", 1.0_dp, aux2, step, 1.0_dp, m, &
1065 0 : filter_eps=eps_filter)
1066 : END IF
1067 :
1068 0 : IF (quadratic_term) THEN
1069 : ! n=step.pq.step
1070 0 : IF (present_oo) THEN
1071 : CALL dbcsr_multiply("N", "N", 1.0_dp, pq, step, 0.0_dp, oo1, &
1072 0 : filter_eps=eps_filter)
1073 : CALL dbcsr_multiply("N", "N", 1.0_dp, oo1, oo, 0.0_dp, oo2, &
1074 0 : filter_eps=eps_filter)
1075 : ELSE
1076 : CALL dbcsr_multiply("N", "N", 1.0_dp, pq, step, 0.0_dp, oo2, &
1077 0 : filter_eps=eps_filter)
1078 : END IF
1079 0 : IF (present_vv) THEN
1080 : CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo2, 0.0_dp, res_trial, &
1081 0 : filter_eps=eps_filter)
1082 : CALL dbcsr_multiply("N", "N", 1.0_dp, vv, res_trial, 0.0_dp, n, &
1083 0 : filter_eps=eps_filter)
1084 : ELSE
1085 : CALL dbcsr_multiply("N", "N", 1.0_dp, step, oo2, 0.0_dp, n, &
1086 0 : filter_eps=eps_filter)
1087 : END IF
1088 :
1089 : ELSE
1090 0 : CALL dbcsr_set(n, 0.0_dp)
1091 : END IF
1092 :
1093 : ! calculate coefficients of the cubic eq for alpha - step size
1094 0 : c0 = 2.0_dp*(dbcsr_frobenius_norm(n))**2
1095 :
1096 0 : CALL dbcsr_dot(m, n, c1)
1097 0 : c1 = -3.0_dp*c1
1098 :
1099 0 : CALL dbcsr_dot(res, n, c2)
1100 0 : c2 = -2.0_dp*c2 + (dbcsr_frobenius_norm(m))**2
1101 :
1102 0 : CALL dbcsr_dot(res, m, c3)
1103 :
1104 : ! find step size
1105 0 : CALL analytic_line_search(c0, c1, c2, c3, step_size, nsteps)
1106 :
1107 0 : IF (nsteps .EQ. 0) THEN
1108 0 : CPABORT("no step sizes!")
1109 : END IF
1110 : ! if we have several possible step sizes
1111 : ! choose one with the lowest objective function
1112 0 : best_norm = 1.0E+100_dp
1113 0 : best_step_size = 0.0_dp
1114 0 : DO istep = 1, nsteps
1115 : ! recompute the residues
1116 0 : CALL dbcsr_copy(res_trial, res)
1117 0 : CALL dbcsr_add(res_trial, m, 1.0_dp, step_size(istep))
1118 0 : CALL dbcsr_add(res_trial, n, 1.0_dp, -step_size(istep)*step_size(istep))
1119 0 : CALL dbcsr_filter(res_trial, eps_filter)
1120 : ! RZK-warning objective function might be different in the case of
1121 : ! tensor_up_down
1122 : !obj_function=0.5_dp*(dbcsr_frobenius_norm(res_trial))**2
1123 0 : obj_function = dbcsr_maxabs(res_trial)
1124 0 : IF (obj_function .LT. best_norm) THEN
1125 0 : best_norm = obj_function
1126 0 : best_step_size = step_size(istep)
1127 : END IF
1128 : END DO
1129 :
1130 : END IF
1131 :
1132 : ! update X along the line
1133 0 : CALL dbcsr_add(x, step, 1.0_dp, best_step_size)
1134 0 : CALL dbcsr_filter(x, eps_filter)
1135 :
1136 : ! evaluate current energy correction
1137 : !change_ecorr=ecorr
1138 : !CALL dbcsr_dot(qp,x,ecorr,"T","N")
1139 : !change_ecorr=ecorr-change_ecorr
1140 :
1141 : ! check convergence and other exit criteria
1142 0 : converged = (best_norm .LT. eps_convergence)
1143 0 : IF (converged .OR. (iteration .GE. max_iter)) THEN
1144 0 : prepare_to_exit = .TRUE.
1145 : END IF
1146 :
1147 0 : t2 = m_walltime()
1148 :
1149 0 : IF (unit_nr > 0) THEN
1150 : WRITE (unit_nr, '(T6,A,1X,I4,1X,E12.3,F8.3)') &
1151 0 : "RICCATI iter ", iteration, best_norm, t2 - t1
1152 : !WRITE(unit_nr,'(T6,A,1X,I4,1X,F15.9,F15.9,E12.3,F8.3)') &
1153 : ! "RICCATI iter ",iteration,ecorr,change_ecorr,best_norm,t2-t1
1154 : END IF
1155 :
1156 0 : t1 = m_walltime()
1157 :
1158 0 : iteration = iteration + 1
1159 :
1160 0 : IF (prepare_to_exit) EXIT
1161 :
1162 : END DO
1163 :
1164 0 : CALL dbcsr_release(aux1)
1165 0 : CALL dbcsr_release(aux2)
1166 0 : CALL dbcsr_release(grad)
1167 0 : CALL dbcsr_release(step)
1168 0 : CALL dbcsr_release(n)
1169 0 : CALL dbcsr_release(m)
1170 0 : CALL dbcsr_release(oo1)
1171 0 : CALL dbcsr_release(oo2)
1172 0 : CALL dbcsr_release(res_trial)
1173 0 : CALL dbcsr_release(vv_step)
1174 0 : CALL dbcsr_release(step_oo)
1175 :
1176 0 : CALL timestop(handle)
1177 :
1178 0 : END SUBROUTINE solve_riccati_equation
1179 :
1180 : ! **************************************************************************************************
1181 : !> \brief Computes a preconditioner from diagonal elements of ~f_oo, ~f_vv
1182 : !> The preconditioner is approximately equal to
1183 : !> prec_ai ~ (e_a - e_i)^(-2)
1184 : !> However, the real expression is more complex
1185 : !> \param prec ...
1186 : !> \param pp ...
1187 : !> \param qq ...
1188 : !> \param eps_filter ...
1189 : !> \par History
1190 : !> 2011.07 created [Rustam Z Khaliullin]
1191 : !> \author Rustam Z Khaliullin
1192 : ! **************************************************************************************************
1193 0 : SUBROUTINE create_preconditioner(prec, pp, qq, eps_filter)
1194 :
1195 : TYPE(dbcsr_type), INTENT(OUT) :: prec
1196 : TYPE(dbcsr_type), INTENT(IN) :: pp, qq
1197 : REAL(KIND=dp), INTENT(IN) :: eps_filter
1198 :
1199 : CHARACTER(len=*), PARAMETER :: routineN = 'create_preconditioner'
1200 :
1201 : INTEGER :: col, handle, hold, iblock_col, &
1202 : iblock_row, mynode, nblkcols_tot, &
1203 : nblkrows_tot, p_nrows, q_nrows, row
1204 : LOGICAL :: tr
1205 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: p_diagonal, q_diagonal
1206 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_new_block
1207 : TYPE(dbcsr_distribution_type) :: dist
1208 : TYPE(dbcsr_type) :: pp_diag, qq_diag, t1, t2, tmp
1209 :
1210 : !LOGICAL, INTENT(IN) :: use_virt_orbs
1211 :
1212 0 : CALL timeset(routineN, handle)
1213 :
1214 : ! ! copy diagonal elements
1215 : ! CALL dbcsr_get_info(pp,nfullrows_total=nrows)
1216 : ! CALL dbcsr_init(pp_diag)
1217 : ! CALL dbcsr_create(pp_diag,template=pp)
1218 : ! ALLOCATE(diagonal(nrows))
1219 : ! CALL dbcsr_get_diag(pp,diagonal)
1220 : ! CALL dbcsr_add_on_diag(pp_diag,1.0_dp)
1221 : ! CALL dbcsr_set_diag(pp_diag,diagonal)
1222 : ! DEALLOCATE(diagonal)
1223 : !
1224 : ! initialize a matrix to 1.0
1225 0 : CALL dbcsr_create(tmp, template=prec)
1226 : ! use an ugly hack to set all elements of tmp to 1
1227 : ! because dbcsr_set does not do it (despite its name)
1228 : !CALL dbcsr_set(tmp,1.0_dp)
1229 0 : CALL dbcsr_get_info(tmp, distribution=dist)
1230 0 : CALL dbcsr_distribution_get(dist, mynode=mynode)
1231 0 : CALL dbcsr_work_create(tmp, work_mutable=.TRUE.)
1232 0 : nblkrows_tot = dbcsr_nblkrows_total(tmp)
1233 0 : nblkcols_tot = dbcsr_nblkcols_total(tmp)
1234 0 : DO row = 1, nblkrows_tot
1235 0 : DO col = 1, nblkcols_tot
1236 0 : tr = .FALSE.
1237 0 : iblock_row = row
1238 0 : iblock_col = col
1239 0 : CALL dbcsr_get_stored_coordinates(tmp, iblock_row, iblock_col, hold)
1240 0 : IF (hold .EQ. mynode) THEN
1241 0 : NULLIFY (p_new_block)
1242 0 : CALL dbcsr_reserve_block2d(tmp, iblock_row, iblock_col, p_new_block)
1243 0 : CPASSERT(ASSOCIATED(p_new_block))
1244 0 : p_new_block(:, :) = 1.0_dp
1245 : END IF ! mynode
1246 : END DO
1247 : END DO
1248 0 : CALL dbcsr_finalize(tmp)
1249 :
1250 : ! copy diagonal elements of pp into cols of a matrix
1251 0 : CALL dbcsr_get_info(pp, nfullrows_total=p_nrows)
1252 0 : CALL dbcsr_create(pp_diag, template=pp)
1253 0 : ALLOCATE (p_diagonal(p_nrows))
1254 0 : CALL dbcsr_get_diag(pp, p_diagonal)
1255 0 : CALL dbcsr_add_on_diag(pp_diag, 1.0_dp)
1256 0 : CALL dbcsr_set_diag(pp_diag, p_diagonal)
1257 : ! RZK-warning is it possible to use dbcsr_scale_by_vector?
1258 : ! or even insert elements directly in the prev cycles
1259 0 : CALL dbcsr_create(t2, template=prec)
1260 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp, pp_diag, &
1261 0 : 0.0_dp, t2, filter_eps=eps_filter)
1262 :
1263 : ! copy diagonal elements qq into rows of a matrix
1264 0 : CALL dbcsr_get_info(qq, nfullrows_total=q_nrows)
1265 0 : CALL dbcsr_create(qq_diag, template=qq)
1266 0 : ALLOCATE (q_diagonal(q_nrows))
1267 0 : CALL dbcsr_get_diag(qq, q_diagonal)
1268 0 : CALL dbcsr_add_on_diag(qq_diag, 1.0_dp)
1269 0 : CALL dbcsr_set_diag(qq_diag, q_diagonal)
1270 0 : CALL dbcsr_set(tmp, 1.0_dp)
1271 0 : CALL dbcsr_create(t1, template=prec)
1272 : CALL dbcsr_multiply("N", "N", 1.0_dp, qq_diag, tmp, &
1273 0 : 0.0_dp, t1, filter_eps=eps_filter)
1274 :
1275 0 : CALL dbcsr_hadamard_product(t1, t2, prec)
1276 0 : CALL dbcsr_release(t1)
1277 0 : CALL dbcsr_scale(prec, 2.0_dp)
1278 :
1279 : ! Get the diagonal of tr(qq).qq
1280 : CALL dbcsr_multiply("T", "N", 1.0_dp, qq, qq, &
1281 : 0.0_dp, qq_diag, retain_sparsity=.TRUE., &
1282 0 : filter_eps=eps_filter)
1283 0 : CALL dbcsr_get_diag(qq_diag, q_diagonal)
1284 0 : CALL dbcsr_set(qq_diag, 0.0_dp)
1285 0 : CALL dbcsr_add_on_diag(qq_diag, 1.0_dp)
1286 0 : CALL dbcsr_set_diag(qq_diag, q_diagonal)
1287 0 : DEALLOCATE (q_diagonal)
1288 0 : CALL dbcsr_set(tmp, 1.0_dp)
1289 : CALL dbcsr_multiply("N", "N", 1.0_dp, qq_diag, tmp, &
1290 0 : 0.0_dp, t2, filter_eps=eps_filter)
1291 0 : CALL dbcsr_release(qq_diag)
1292 0 : CALL dbcsr_add(prec, t2, 1.0_dp, 1.0_dp)
1293 :
1294 : ! Get the diagonal of pp.tr(pp)
1295 : CALL dbcsr_multiply("N", "T", 1.0_dp, pp, pp, &
1296 : 0.0_dp, pp_diag, retain_sparsity=.TRUE., &
1297 0 : filter_eps=eps_filter)
1298 0 : CALL dbcsr_get_diag(pp_diag, p_diagonal)
1299 0 : CALL dbcsr_set(pp_diag, 0.0_dp)
1300 0 : CALL dbcsr_add_on_diag(pp_diag, 1.0_dp)
1301 0 : CALL dbcsr_set_diag(pp_diag, p_diagonal)
1302 0 : DEALLOCATE (p_diagonal)
1303 0 : CALL dbcsr_set(tmp, 1.0_dp)
1304 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp, pp_diag, &
1305 0 : 0.0_dp, t2, filter_eps=eps_filter)
1306 0 : CALL dbcsr_release(tmp)
1307 0 : CALL dbcsr_release(pp_diag)
1308 0 : CALL dbcsr_add(prec, t2, 1.0_dp, 1.0_dp)
1309 :
1310 : ! now add the residual component
1311 : !CALL dbcsr_hadamard_product(res,qp,t2)
1312 : !CALL dbcsr_add(prec,t2,1.0_dp,-2.0_dp)
1313 0 : CALL dbcsr_release(t2)
1314 0 : CALL inverse_of_elements(prec)
1315 0 : CALL dbcsr_filter(prec, eps_filter)
1316 :
1317 0 : CALL timestop(handle)
1318 :
1319 0 : END SUBROUTINE create_preconditioner
1320 :
1321 : ! **************************************************************************************************
1322 : !> \brief Computes 1/x of the matrix elements.
1323 : !> \param matrix ...
1324 : !> \author Ole Schuett
1325 : ! **************************************************************************************************
1326 0 : SUBROUTINE inverse_of_elements(matrix)
1327 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1328 :
1329 : CHARACTER(len=*), PARAMETER :: routineN = 'inverse_of_elements'
1330 :
1331 : INTEGER :: handle
1332 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: block
1333 : TYPE(dbcsr_iterator_type) :: iter
1334 :
1335 0 : CALL timeset(routineN, handle)
1336 0 : CALL dbcsr_iterator_start(iter, matrix)
1337 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1338 0 : CALL dbcsr_iterator_next_block(iter, block=block)
1339 0 : block = 1.0_dp/block
1340 : END DO
1341 0 : CALL dbcsr_iterator_stop(iter)
1342 0 : CALL timestop(handle)
1343 :
1344 0 : END SUBROUTINE inverse_of_elements
1345 :
1346 : ! **************************************************************************************************
1347 : !> \brief Finds real roots of a cubic equation
1348 : !> > a*x**3 + b*x**2 + c*x + d = 0
1349 : !> and returns only those roots for which the derivative is positive
1350 : !>
1351 : !> Step 0: Check the true order of the equation. Cubic, quadratic, linear?
1352 : !> Step 1: Calculate p and q
1353 : !> p = ( 3*c/a - (b/a)**2 ) / 3
1354 : !> q = ( 2*(b/a)**3 - 9*b*c/a/a + 27*d/a ) / 27
1355 : !> Step 2: Calculate discriminant D
1356 : !> D = (p/3)**3 + (q/2)**2
1357 : !> Step 3: Depending on the sign of D, we follow different strategy.
1358 : !> If D<0, three distinct real roots.
1359 : !> If D=0, three real roots of which at least two are equal.
1360 : !> If D>0, one real and two complex roots.
1361 : !> Step 3a: For D>0 and D=0,
1362 : !> Calculate u and v
1363 : !> u = cubic_root(-q/2 + sqrt(D))
1364 : !> v = cubic_root(-q/2 - sqrt(D))
1365 : !> Find the three transformed roots
1366 : !> y1 = u + v
1367 : !> y2 = -(u+v)/2 + i (u-v)*sqrt(3)/2
1368 : !> y3 = -(u+v)/2 - i (u-v)*sqrt(3)/2
1369 : !> Step 3b Alternately, for D<0, a trigonometric formulation is more convenient
1370 : !> y1 = 2 * sqrt(|p|/3) * cos(phi/3)
1371 : !> y2 = -2 * sqrt(|p|/3) * cos((phi+pi)/3)
1372 : !> y3 = -2 * sqrt(|p|/3) * cos((phi-pi)/3)
1373 : !> where phi = acos(-q/2/sqrt(|p|**3/27))
1374 : !> pi = 3.141592654...
1375 : !> Step 4 Find the real roots
1376 : !> x = y - b/a/3
1377 : !> Step 5 Check the derivative and return only those real roots
1378 : !> for which the derivative is positive
1379 : !>
1380 : !> \param a ...
1381 : !> \param b ...
1382 : !> \param c ...
1383 : !> \param d ...
1384 : !> \param minima ...
1385 : !> \param nmins ...
1386 : !> \par History
1387 : !> 2011.06 created [Rustam Z Khaliullin]
1388 : !> \author Rustam Z Khaliullin
1389 : ! **************************************************************************************************
1390 0 : SUBROUTINE analytic_line_search(a, b, c, d, minima, nmins)
1391 :
1392 : REAL(KIND=dp), INTENT(IN) :: a, b, c, d
1393 : REAL(KIND=dp), DIMENSION(3), INTENT(OUT) :: minima
1394 : INTEGER, INTENT(OUT) :: nmins
1395 :
1396 : INTEGER :: i, nroots
1397 : REAL(KIND=dp) :: DD, der, p, phi, pi, q, temp1, temp2, u, &
1398 : v, y1, y2, y2i, y2r, y3
1399 : REAL(KIND=dp), DIMENSION(3) :: x
1400 :
1401 : ! CALL timeset(routineN,handle)
1402 :
1403 0 : pi = ACOS(-1.0_dp)
1404 :
1405 : ! Step 0: Check coefficients and find the true order of the eq
1406 0 : IF (a .EQ. 0.0_dp) THEN
1407 0 : IF (b .EQ. 0.0_dp) THEN
1408 0 : IF (c .EQ. 0.0_dp) THEN
1409 : ! Non-equation, no valid solutions
1410 : nroots = 0
1411 : ELSE
1412 : ! Linear equation with one root.
1413 0 : nroots = 1
1414 0 : x(1) = -d/c
1415 : END IF
1416 : ELSE
1417 : ! Quadratic equation with max two roots.
1418 0 : DD = c*c - 4.0_dp*b*d
1419 0 : IF (DD .GT. 0.0_dp) THEN
1420 0 : nroots = 2
1421 0 : x(1) = (-c + SQRT(DD))/2.0_dp/b
1422 0 : x(2) = (-c - SQRT(DD))/2.0_dp/b
1423 0 : ELSE IF (DD .LT. 0.0_dp) THEN
1424 : nroots = 0
1425 : ELSE
1426 0 : nroots = 1
1427 0 : x(1) = -c/2.0_dp/b
1428 : END IF
1429 : END IF
1430 : ELSE
1431 : ! Cubic equation with max three roots
1432 : ! Calculate p and q
1433 0 : p = c/a - b*b/a/a/3.0_dp
1434 0 : q = (2.0_dp*b*b*b/a/a/a - 9.0_dp*b*c/a/a + 27.0_dp*d/a)/27.0_dp
1435 :
1436 : ! Calculate DD
1437 0 : DD = p*p*p/27.0_dp + q*q/4.0_dp
1438 :
1439 0 : IF (DD .LT. 0.0_dp) THEN
1440 : ! three real unequal roots -- use the trigonometric formulation
1441 0 : phi = ACOS(-q/2.0_dp/SQRT(ABS(p*p*p)/27.0_dp))
1442 0 : temp1 = 2.0_dp*SQRT(ABS(p)/3.0_dp)
1443 0 : y1 = temp1*COS(phi/3.0_dp)
1444 0 : y2 = -temp1*COS((phi + pi)/3.0_dp)
1445 0 : y3 = -temp1*COS((phi - pi)/3.0_dp)
1446 : ELSE
1447 : ! 1 real & 2 conjugate complex roots OR 3 real roots (some are equal)
1448 0 : temp1 = -q/2.0_dp + SQRT(DD)
1449 0 : temp2 = -q/2.0_dp - SQRT(DD)
1450 0 : u = ABS(temp1)**(1.0_dp/3.0_dp)
1451 0 : v = ABS(temp2)**(1.0_dp/3.0_dp)
1452 0 : IF (temp1 .LT. 0.0_dp) u = -u
1453 0 : IF (temp2 .LT. 0.0_dp) v = -v
1454 0 : y1 = u + v
1455 0 : y2r = -(u + v)/2.0_dp
1456 0 : y2i = (u - v)*SQRT(3.0_dp)/2.0_dp
1457 : END IF
1458 :
1459 : ! Final transformation
1460 0 : temp1 = b/a/3.0_dp
1461 0 : y1 = y1 - temp1
1462 0 : y2 = y2 - temp1
1463 0 : y3 = y3 - temp1
1464 0 : y2r = y2r - temp1
1465 :
1466 : ! Assign answers
1467 0 : IF (DD .LT. 0.0_dp) THEN
1468 0 : nroots = 3
1469 0 : x(1) = y1
1470 0 : x(2) = y2
1471 0 : x(3) = y3
1472 0 : ELSE IF (DD .EQ. 0.0_dp) THEN
1473 0 : nroots = 2
1474 0 : x(1) = y1
1475 0 : x(2) = y2r
1476 : !x(3) = cmplx(y2r, 0.)
1477 : ELSE
1478 0 : nroots = 1
1479 0 : x(1) = y1
1480 : !x(2) = cmplx(y2r, y2i)
1481 : !x(3) = cmplx(y2r,-y2i)
1482 : END IF
1483 :
1484 : END IF
1485 :
1486 : !write(*,'(i2,a)') nroots, ' real root(s)'
1487 0 : nmins = 0
1488 0 : DO i = 1, nroots
1489 : ! maximum or minimum? use the derivative
1490 : ! 3*a*x**2+2*b*x+c
1491 0 : der = 3.0_dp*a*x(i)*x(i) + 2.0_dp*b*x(i) + c
1492 0 : IF (der .GT. 0.0_dp) THEN
1493 0 : nmins = nmins + 1
1494 0 : minima(nmins) = x(i)
1495 : !write(*,'(a,i2,a,f10.5)') 'Minimum ', i, ', value: ', x(i)
1496 : END IF
1497 : END DO
1498 :
1499 : ! CALL timestop(handle)
1500 :
1501 0 : END SUBROUTINE analytic_line_search
1502 :
1503 : ! **************************************************************************************************
1504 : !> \brief Diagonalizes diagonal blocks of a symmetric dbcsr matrix
1505 : !> and returs its eigenvectors
1506 : !> \param matrix ...
1507 : !> \param c ...
1508 : !> \param e ...
1509 : !> \par History
1510 : !> 2011.07 created [Rustam Z Khaliullin]
1511 : !> \author Rustam Z Khaliullin
1512 : ! **************************************************************************************************
1513 0 : SUBROUTINE diagonalize_diagonal_blocks(matrix, c, e)
1514 :
1515 : TYPE(dbcsr_type), INTENT(IN) :: matrix
1516 : TYPE(dbcsr_type), INTENT(OUT) :: c
1517 : TYPE(dbcsr_type), INTENT(OUT), OPTIONAL :: e
1518 :
1519 : CHARACTER(len=*), PARAMETER :: routineN = 'diagonalize_diagonal_blocks'
1520 :
1521 : INTEGER :: handle, iblock_col, iblock_row, &
1522 : iblock_size, info, lwork, orbital
1523 : LOGICAL :: block_needed, do_eigenvalues
1524 0 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: eigenvalues, work
1525 0 : REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: data_copy
1526 0 : REAL(kind=dp), DIMENSION(:, :), POINTER :: data_p, p_new_block
1527 : TYPE(dbcsr_iterator_type) :: iter
1528 :
1529 0 : CALL timeset(routineN, handle)
1530 :
1531 0 : IF (PRESENT(e)) THEN
1532 : do_eigenvalues = .TRUE.
1533 : ELSE
1534 0 : do_eigenvalues = .FALSE.
1535 : END IF
1536 :
1537 : ! create a matrix for eigenvectors
1538 0 : CALL dbcsr_work_create(c, work_mutable=.TRUE.)
1539 0 : IF (do_eigenvalues) &
1540 0 : CALL dbcsr_work_create(e, work_mutable=.TRUE.)
1541 :
1542 0 : CALL dbcsr_iterator_start(iter, matrix)
1543 :
1544 0 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1545 :
1546 0 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, data_p, row_size=iblock_size)
1547 :
1548 0 : block_needed = .FALSE.
1549 0 : IF (iblock_row == iblock_col) block_needed = .TRUE.
1550 :
1551 0 : IF (block_needed) THEN
1552 :
1553 : ! Prepare data
1554 0 : ALLOCATE (eigenvalues(iblock_size))
1555 0 : ALLOCATE (data_copy(iblock_size, iblock_size))
1556 0 : data_copy(:, :) = data_p(:, :)
1557 :
1558 : ! Query the optimal workspace for dsyev
1559 0 : LWORK = -1
1560 0 : ALLOCATE (WORK(MAX(1, LWORK)))
1561 0 : CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
1562 0 : LWORK = INT(WORK(1))
1563 0 : DEALLOCATE (WORK)
1564 :
1565 : ! Allocate the workspace and solve the eigenproblem
1566 0 : ALLOCATE (WORK(MAX(1, LWORK)))
1567 0 : CALL dsyev('V', 'L', iblock_size, data_copy, iblock_size, eigenvalues, WORK, LWORK, INFO)
1568 0 : IF (INFO .NE. 0) CPABORT("DSYEV failed")
1569 :
1570 : ! copy eigenvectors into a cp_dbcsr matrix
1571 0 : NULLIFY (p_new_block)
1572 0 : CALL dbcsr_reserve_block2d(c, iblock_row, iblock_col, p_new_block)
1573 0 : CPASSERT(ASSOCIATED(p_new_block))
1574 0 : p_new_block(:, :) = data_copy(:, :)
1575 :
1576 : ! if requested copy eigenvalues into a cp_dbcsr matrix
1577 0 : IF (do_eigenvalues) THEN
1578 0 : NULLIFY (p_new_block)
1579 0 : CALL dbcsr_reserve_block2d(e, iblock_row, iblock_col, p_new_block)
1580 0 : CPASSERT(ASSOCIATED(p_new_block))
1581 0 : p_new_block(:, :) = 0.0_dp
1582 0 : DO orbital = 1, iblock_size
1583 0 : p_new_block(orbital, orbital) = eigenvalues(orbital)
1584 : END DO
1585 : END IF
1586 :
1587 0 : DEALLOCATE (WORK)
1588 0 : DEALLOCATE (data_copy)
1589 0 : DEALLOCATE (eigenvalues)
1590 :
1591 : END IF
1592 :
1593 : END DO
1594 :
1595 0 : CALL dbcsr_iterator_stop(iter)
1596 :
1597 0 : CALL dbcsr_finalize(c)
1598 0 : IF (do_eigenvalues) CALL dbcsr_finalize(e)
1599 :
1600 0 : CALL timestop(handle)
1601 :
1602 0 : END SUBROUTINE diagonalize_diagonal_blocks
1603 :
1604 : ! **************************************************************************************************
1605 : !> \brief Transforms a matrix M_out = tr(U1) * M_in * U2
1606 : !> \param matrix ...
1607 : !> \param u1 ...
1608 : !> \param u2 ...
1609 : !> \param eps_filter ...
1610 : !> \par History
1611 : !> 2011.10 created [Rustam Z Khaliullin]
1612 : !> \author Rustam Z Khaliullin
1613 : ! **************************************************************************************************
1614 0 : SUBROUTINE matrix_forward_transform(matrix, u1, u2, eps_filter)
1615 :
1616 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1617 : TYPE(dbcsr_type), INTENT(IN) :: u1, u2
1618 : REAL(KIND=dp), INTENT(IN) :: eps_filter
1619 :
1620 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_forward_transform'
1621 :
1622 : INTEGER :: handle
1623 : TYPE(dbcsr_type) :: tmp
1624 :
1625 0 : CALL timeset(routineN, handle)
1626 :
1627 : CALL dbcsr_create(tmp, template=matrix, &
1628 0 : matrix_type=dbcsr_type_no_symmetry)
1629 : CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, u2, 0.0_dp, tmp, &
1630 0 : filter_eps=eps_filter)
1631 : CALL dbcsr_multiply("T", "N", 1.0_dp, u1, tmp, 0.0_dp, matrix, &
1632 0 : filter_eps=eps_filter)
1633 0 : CALL dbcsr_release(tmp)
1634 :
1635 0 : CALL timestop(handle)
1636 :
1637 0 : END SUBROUTINE matrix_forward_transform
1638 :
1639 : ! **************************************************************************************************
1640 : !> \brief Transforms a matrix M_out = U1 * M_in * tr(U2)
1641 : !> \param matrix ...
1642 : !> \param u1 ...
1643 : !> \param u2 ...
1644 : !> \param eps_filter ...
1645 : !> \par History
1646 : !> 2011.10 created [Rustam Z Khaliullin]
1647 : !> \author Rustam Z Khaliullin
1648 : ! **************************************************************************************************
1649 0 : SUBROUTINE matrix_backward_transform(matrix, u1, u2, eps_filter)
1650 :
1651 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix
1652 : TYPE(dbcsr_type), INTENT(IN) :: u1, u2
1653 : REAL(KIND=dp), INTENT(IN) :: eps_filter
1654 :
1655 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_backward_transform'
1656 :
1657 : INTEGER :: handle
1658 : TYPE(dbcsr_type) :: tmp
1659 :
1660 0 : CALL timeset(routineN, handle)
1661 :
1662 : CALL dbcsr_create(tmp, template=matrix, &
1663 0 : matrix_type=dbcsr_type_no_symmetry)
1664 : CALL dbcsr_multiply("N", "T", 1.0_dp, matrix, u2, 0.0_dp, tmp, &
1665 0 : filter_eps=eps_filter)
1666 : CALL dbcsr_multiply("N", "N", 1.0_dp, u1, tmp, 0.0_dp, matrix, &
1667 0 : filter_eps=eps_filter)
1668 0 : CALL dbcsr_release(tmp)
1669 :
1670 0 : CALL timestop(handle)
1671 :
1672 0 : END SUBROUTINE matrix_backward_transform
1673 :
1674 : !! **************************************************************************************************
1675 : !!> \brief Transforms to a representation in which diagonal blocks
1676 : !!> of qq and pp matrices are diagonal. This can improve convergence
1677 : !!> of PCG
1678 : !!> \par History
1679 : !!> 2011.07 created [Rustam Z Khaliullin]
1680 : !!> \author Rustam Z Khaliullin
1681 : !! **************************************************************************************************
1682 : ! SUBROUTINE transform_matrices_to_blk_diag(matrix_pp,matrix_qq,matrix_qp,&
1683 : ! matrix_pq,eps_filter)
1684 : !
1685 : ! TYPE(dbcsr_type), INTENT(INOUT) :: matrix_pp, matrix_qq,&
1686 : ! matrix_qp, matrix_pq
1687 : ! REAL(KIND=dp), INTENT(IN) :: eps_filter
1688 : !
1689 : ! CHARACTER(len=*), PARAMETER :: routineN = 'transform_matrices_to_blk_diag',&
1690 : ! routineP = moduleN//':'//routineN
1691 : !
1692 : ! TYPE(dbcsr_type) :: tmp_pp, tmp_qq,&
1693 : ! tmp_qp, tmp_pq,&
1694 : ! blk, blk2
1695 : ! INTEGER :: handle
1696 : !
1697 : ! CALL timeset(routineN,handle)
1698 : !
1699 : ! ! find a better basis by diagonalizing diagonal blocks
1700 : ! ! first pp
1701 : ! CALL dbcsr_init(blk)
1702 : ! CALL dbcsr_create(blk,template=matrix_pp)
1703 : ! CALL diagonalize_diagonal_blocks(matrix_pp,blk)
1704 : !
1705 : ! ! convert matrices to the new basis
1706 : ! CALL dbcsr_init(tmp_pp)
1707 : ! CALL dbcsr_create(tmp_pp,template=matrix_pp)
1708 : ! CALL dbcsr_multiply("N","N",1.0_dp,matrix_pp,blk,0.0_dp,tmp_pp,&
1709 : ! filter_eps=eps_filter)
1710 : ! CALL dbcsr_multiply("T","N",1.0_dp,blk,tmp_pp,0.0_dp,matrix_pp,&
1711 : ! filter_eps=eps_filter)
1712 : ! CALL dbcsr_release(tmp_pp)
1713 : !
1714 : ! ! now qq
1715 : ! CALL dbcsr_init(blk2)
1716 : ! CALL dbcsr_create(blk2,template=matrix_qq)
1717 : ! CALL diagonalize_diagonal_blocks(matrix_qq,blk2)
1718 : !
1719 : ! CALL dbcsr_init(tmp_qq)
1720 : ! CALL dbcsr_create(tmp_qq,template=matrix_qq)
1721 : ! CALL dbcsr_multiply("N","N",1.0_dp,matrix_qq,blk2,0.0_dp,tmp_qq,&
1722 : ! filter_eps=eps_filter)
1723 : ! CALL dbcsr_multiply("T","N",1.0_dp,blk2,tmp_qq,0.0_dp,matrix_qq,&
1724 : ! filter_eps=eps_filter)
1725 : ! CALL dbcsr_release(tmp_qq)
1726 : !
1727 : ! ! transform pq
1728 : ! CALL dbcsr_init(tmp_pq)
1729 : ! CALL dbcsr_create(tmp_pq,template=matrix_pq)
1730 : ! CALL dbcsr_multiply("T","N",1.0_dp,blk,matrix_pq,0.0_dp,tmp_pq,&
1731 : ! filter_eps=eps_filter)
1732 : ! CALL dbcsr_multiply("N","N",1.0_dp,tmp_pq,blk2,0.0_dp,matrix_pq,&
1733 : ! filter_eps=eps_filter)
1734 : ! CALL dbcsr_release(tmp_pq)
1735 : !
1736 : ! ! transform qp
1737 : ! CALL dbcsr_init(tmp_qp)
1738 : ! CALL dbcsr_create(tmp_qp,template=matrix_qp)
1739 : ! CALL dbcsr_multiply("N","N",1.0_dp,matrix_qp,blk,0.0_dp,tmp_qp,&
1740 : ! filter_eps=eps_filter)
1741 : ! CALL dbcsr_multiply("T","N",1.0_dp,blk2,tmp_qp,0.0_dp,matrix_qp,&
1742 : ! filter_eps=eps_filter)
1743 : ! CALL dbcsr_release(tmp_qp)
1744 : !
1745 : ! CALL dbcsr_release(blk2)
1746 : ! CALL dbcsr_release(blk)
1747 : !
1748 : ! CALL timestop(handle)
1749 : !
1750 : ! END SUBROUTINE transform_matrices_to_blk_diag
1751 :
1752 : ! **************************************************************************************************
1753 : !> \brief computes oo, ov, vo, and vv blocks of the ks matrix
1754 : !> \par History
1755 : !> 2011.06 created [Rustam Z Khaliullin]
1756 : !> \author Rustam Z Khaliullin
1757 : ! **************************************************************************************************
1758 : ! SUBROUTINE ct_step_env_execute(env)
1759 : !
1760 : ! TYPE(ct_step_env_type) :: env
1761 : !
1762 : ! CHARACTER(len=*), PARAMETER :: routineN = 'ct_step_env_execute', &
1763 : ! routineP = moduleN//':'//routineN
1764 : !
1765 : ! INTEGER :: handle
1766 : !
1767 : ! CALL timeset(routineN,handle)
1768 : !
1769 : !
1770 : ! CALL timestop(handle)
1771 : !
1772 : ! END SUBROUTINE ct_step_env_execute
1773 :
1774 : END MODULE ct_methods
1775 :
|
