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