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