(git:ba1d7ca)
Loading...
Searching...
No Matches
ec_orth_solver.F
Go to the documentation of this file.
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 AO-based conjugate-gradient response solver routines
10!>
11!>
12!> \date 09.2019
13!> \author Fabian Belleflamme
14! **************************************************************************************************
15MODULE ec_orth_solver
16 USE admm_types, ONLY: admm_type,&
17 get_admm_env
18 USE cp_control_types, ONLY: dft_control_type
19 USE cp_dbcsr_api, ONLY: &
20 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_filter, dbcsr_finalize, &
21 dbcsr_get_info, dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_scale, dbcsr_set, &
22 dbcsr_transposed, dbcsr_type, dbcsr_type_no_symmetry
23 USE cp_dbcsr_contrib, ONLY: dbcsr_add_on_diag,&
24 dbcsr_checksum,&
25 dbcsr_dot
26 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
27 dbcsr_deallocate_matrix_set
28 USE cp_external_control, ONLY: external_control
29 USE input_constants, ONLY: do_admm_aux_exch_func_none,&
30 kg_tnadd_embed,&
31 kg_tnadd_embed_ri,&
32 ls_s_sqrt_ns,&
33 ls_s_sqrt_proot,&
34 precond_mlp
35 USE input_section_types, ONLY: section_vals_get_subs_vals,&
36 section_vals_type,&
37 section_vals_val_get
38 USE iterate_matrix, ONLY: matrix_sqrt_newton_schulz,&
39 matrix_sqrt_proot
40 USE kg_correction, ONLY: kg_ekin_subset
41 USE kinds, ONLY: dp
42 USE machine, ONLY: m_flush,&
43 m_walltime
44 USE mathlib, ONLY: abnormal_value
45 USE message_passing, ONLY: mp_para_env_type
46 USE pw_env_types, ONLY: pw_env_get,&
47 pw_env_type
48 USE pw_methods, ONLY: pw_axpy,&
49 pw_scale,&
50 pw_transfer,&
51 pw_zero
52 USE pw_poisson_methods, ONLY: pw_poisson_solve
53 USE pw_poisson_types, ONLY: pw_poisson_type
54 USE pw_pool_types, ONLY: pw_pool_p_type,&
55 pw_pool_type
56 USE pw_types, ONLY: pw_c1d_gs_type,&
57 pw_r3d_rs_type
58 USE qs_environment_types, ONLY: get_qs_env,&
59 qs_environment_type
60 USE qs_fxc, ONLY: qs_fxc_create,&
61 qs_fxc_prep
62 USE qs_integrate_potential, ONLY: integrate_v_rspace
63 USE qs_kpp1_env_types, ONLY: qs_kpp1_env_type
64 USE qs_linres_kernel, ONLY: apply_hfx,&
65 apply_xc_admm
66 USE qs_linres_types, ONLY: linres_control_type
67 USE qs_p_env_methods, ONLY: p_env_check_i_alloc,&
68 p_env_finish_kpp1,&
69 p_env_update_rho
70 USE qs_p_env_types, ONLY: qs_p_env_type
71 USE qs_rho_atom_types, ONLY: rho_atom_type
72 USE qs_rho_types, ONLY: qs_rho_get,&
73 qs_rho_type
74#include "./base/base_uses.f90"
75
76 IMPLICIT NONE
77
78 PRIVATE
79
80 ! Global parameters
81
82 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ec_orth_solver'
83
84 ! Public subroutines
85
86 PUBLIC :: ec_response_ao
87
88CONTAINS
89
90! **************************************************************************************************
91!> \brief Preconditioning of the AO-based CG linear response solver
92!> M * z_0 = r_0
93!> M(X) = [F,B], with B = [X,P]
94!> for M we need F and P in ortho basis
95!> Returns z_0, the preconditioned residual in orthonormal basis
96!>
97!> All matrices are in orthonormal Lowdin basis
98!>
99!> \param qs_env ...
100!> \param matrix_ks Ground-state Kohn-Sham matrix
101!> \param matrix_p Ground-state Density matrix
102!> \param matrix_rhs Unpreconditioned residual of linear response CG
103!> \param matrix_cg_z Preconditioned residual
104!> \param eps_filter ...
105!> \param iounit ...
106!>
107!> \param silent ...
108!> \date 01.2020
109!> \author Fabian Belleflamme
110! **************************************************************************************************
111 SUBROUTINE ec_preconditioner(qs_env, matrix_ks, matrix_p, matrix_rhs, &
112 matrix_cg_z, eps_filter, iounit, silent)
113
114 TYPE(qs_environment_type), POINTER :: qs_env
115 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
116 POINTER :: matrix_ks, matrix_p, matrix_rhs
117 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
118 POINTER :: matrix_cg_z
119 REAL(KIND=dp), INTENT(IN) :: eps_filter
120 INTEGER, INTENT(IN) :: iounit
121 LOGICAL, INTENT(IN), OPTIONAL :: silent
122
123 CHARACTER(len=*), PARAMETER :: routineN = 'ec_preconditioner'
124
125 INTEGER :: handle, i, ispin, max_iter, nao, nspins
126 LOGICAL :: converged, my_silent
127 REAL(KIND=dp) :: norm_res, t1, t2
128 REAL(KIND=dp), DIMENSION(:), POINTER :: alpha, beta, new_norm, norm_ca, norm_rr
129 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_Ax, matrix_b, matrix_cg, &
130 matrix_res
131 TYPE(dft_control_type), POINTER :: dft_control
132 TYPE(linres_control_type), POINTER :: linres_control
133
134 CALL timeset(routinen, handle)
135
136 my_silent = .false.
137 IF (PRESENT(silent)) my_silent = silent
138
139 cpassert(ASSOCIATED(qs_env))
140 cpassert(ASSOCIATED(matrix_ks))
141 cpassert(ASSOCIATED(matrix_p))
142 cpassert(ASSOCIATED(matrix_rhs))
143 cpassert(ASSOCIATED(matrix_cg_z))
144
145 NULLIFY (dft_control, linres_control)
146
147 t1 = m_walltime()
148
149 CALL get_qs_env(qs_env=qs_env, &
150 dft_control=dft_control, &
151 linres_control=linres_control)
152 nspins = dft_control%nspins
153 CALL dbcsr_get_info(matrix_ks(1)%matrix, nfullrows_total=nao)
154
155 ALLOCATE (alpha(nspins), beta(nspins), new_norm(nspins), norm_ca(nspins), norm_rr(nspins))
156
157 !----------------------------------------
158 ! Create non-symmetric matrices: Ax, B, cg, res
159 !----------------------------------------
160
161 NULLIFY (matrix_ax, matrix_b, matrix_cg, matrix_res)
162 CALL dbcsr_allocate_matrix_set(matrix_ax, nspins)
163 CALL dbcsr_allocate_matrix_set(matrix_b, nspins)
164 CALL dbcsr_allocate_matrix_set(matrix_cg, nspins)
165 CALL dbcsr_allocate_matrix_set(matrix_res, nspins)
166
167 DO ispin = 1, nspins
168 ALLOCATE (matrix_ax(ispin)%matrix)
169 ALLOCATE (matrix_b(ispin)%matrix)
170 ALLOCATE (matrix_cg(ispin)%matrix)
171 ALLOCATE (matrix_res(ispin)%matrix)
172 CALL dbcsr_create(matrix_ax(ispin)%matrix, name="linop MATRIX", &
173 template=matrix_ks(1)%matrix, &
174 matrix_type=dbcsr_type_no_symmetry)
175 CALL dbcsr_create(matrix_b(ispin)%matrix, name="MATRIX B", &
176 template=matrix_ks(1)%matrix, &
177 matrix_type=dbcsr_type_no_symmetry)
178 CALL dbcsr_create(matrix_cg(ispin)%matrix, name="TRIAL MATRIX", &
179 template=matrix_ks(1)%matrix, &
180 matrix_type=dbcsr_type_no_symmetry)
181 CALL dbcsr_create(matrix_res(ispin)%matrix, name="RESIDUE", &
182 template=matrix_ks(1)%matrix, &
183 matrix_type=dbcsr_type_no_symmetry)
184 END DO
185
186 !----------------------------------------
187 ! Get righ-hand-side operators
188 !----------------------------------------
189
190 ! Initial guess z_0
191 DO ispin = 1, nspins
192 CALL dbcsr_copy(matrix_cg_z(ispin)%matrix, matrix_rhs(ispin)%matrix)
193
194 ! r_0 = b
195 CALL dbcsr_copy(matrix_res(ispin)%matrix, matrix_rhs(ispin)%matrix)
196 END DO
197
198 ! Projector on trial matrix
199 ! Projector does not need to be applied here,
200 ! as matrix_rhs already had this done before entering preconditioner
201 !CALL projector(qs_env, matrix_p, matrix_cg_z, eps_filter)
202
203 ! Mz_0
204 CALL hessian_op1(matrix_ks, matrix_p, matrix_cg_z, matrix_b, matrix_ax, eps_filter)
205
206 ! r_0 = b - Ax_0
207 DO ispin = 1, nspins
208 CALL dbcsr_add(matrix_res(ispin)%matrix, matrix_ax(ispin)%matrix, 1.0_dp, -1.0_dp)
209 END DO
210
211 ! Matrix projector T
212 CALL projector(qs_env, matrix_p, matrix_res, eps_filter)
213
214 DO ispin = 1, nspins
215 ! cg = p_0 = z_0
216 CALL dbcsr_copy(matrix_cg(ispin)%matrix, matrix_res(ispin)%matrix)
217 END DO
218
219 ! header
220 IF (iounit > 0 .AND. .NOT. my_silent) THEN
221 WRITE (iounit, "(/,T10,A)") "Preconditioning of search direction"
222 WRITE (iounit, "(/,T10,A,T25,A,T42,A,T62,A,/,T10,A)") &
223 "Iteration", "Stepsize", "Convergence", "Time", &
224 repeat("-", 58)
225 END IF
226
227 alpha(:) = 0.0_dp
228 max_iter = 200
229 converged = .false.
230 norm_res = 0.0_dp
231
232 ! start iteration
233 iteration: DO i = 1, max_iter
234
235 ! Hessian Ax = [F,B] is updated preconditioner
236 CALL hessian_op1(matrix_ks, matrix_p, matrix_cg, matrix_b, matrix_ax, eps_filter)
237
238 ! Matrix projector
239 CALL projector(qs_env, matrix_p, matrix_ax, eps_filter)
240
241 DO ispin = 1, nspins
242
243 ! Tr(r_0 * r_0)
244 CALL dbcsr_dot(matrix_res(ispin)%matrix, matrix_res(ispin)%matrix, norm_rr(ispin))
245 IF (abnormal_value(norm_rr(ispin))) THEN
246 cpabort("Preconditioner: Tr[r_j*r_j] is an abnormal value (NaN/Inf)")
247 END IF
248
249 IF (norm_rr(ispin) < 0.0_dp) cpabort("norm_rr < 0")
250 norm_res = max(norm_res, abs(norm_rr(ispin)/real(nao, dp)))
251
252 ! norm_cA = tr(Ap_j * p_j)
253 CALL dbcsr_dot(matrix_cg(ispin)%matrix, matrix_ax(ispin)%matrix, norm_ca(ispin))
254
255 ! Determine step-size
256 IF (norm_ca(ispin) < linres_control%eps) THEN
257 alpha(ispin) = 1.0_dp
258 ELSE
259 alpha(ispin) = norm_rr(ispin)/norm_ca(ispin)
260 END IF
261
262 ! x_j+1 = x_j + alpha*p_j
263 ! save contribution of this iteration
264 CALL dbcsr_add(matrix_cg_z(ispin)%matrix, matrix_cg(ispin)%matrix, 1.0_dp, alpha(ispin))
265
266 ! r_j+1 = r_j - alpha * Ap_j
267 CALL dbcsr_add(matrix_res(ispin)%matrix, matrix_ax(ispin)%matrix, 1.0_dp, -alpha(ispin))
268
269 END DO
270
271 norm_res = 0.0_dp
272
273 DO ispin = 1, nspins
274 ! Tr[r_j+1*z_j+1]
275 CALL dbcsr_dot(matrix_res(ispin)%matrix, matrix_res(ispin)%matrix, new_norm(ispin))
276 IF (new_norm(ispin) < 0.0_dp) cpabort("tr(r_j+1*z_j+1) < 0")
277 IF (abnormal_value(new_norm(ispin))) THEN
278 cpabort("Preconditioner: Tr[r_j+1*z_j+1] is an abnormal value (NaN/Inf)")
279 END IF
280 norm_res = max(norm_res, new_norm(ispin)/real(nao, dp))
281
282 IF (norm_rr(ispin) < linres_control%eps*0.001_dp &
283 .OR. new_norm(ispin) < linres_control%eps*0.001_dp) THEN
284 beta(ispin) = 0.0_dp
285 converged = .true.
286 ELSE
287 beta(ispin) = new_norm(ispin)/norm_rr(ispin)
288 END IF
289
290 ! update new search vector (matrix cg)
291 ! cg_j+1 = z_j+1 + beta*cg_j
292 CALL dbcsr_add(matrix_cg(ispin)%matrix, matrix_res(ispin)%matrix, beta(ispin), 1.0_dp)
293 CALL dbcsr_filter(matrix_cg(ispin)%matrix, eps_filter)
294
295 norm_rr(ispin) = new_norm(ispin)
296 END DO
297
298 ! Convergence criteria
299 IF (norm_res < linres_control%eps) THEN
300 converged = .true.
301 END IF
302
303 t2 = m_walltime()
304 IF (i == 1 .OR. mod(i, 1) == 0 .OR. converged) THEN
305 IF (iounit > 0 .AND. .NOT. my_silent) THEN
306 WRITE (iounit, "(T10,I5,T25,1E8.2,T33,F25.14,T58,F8.2)") &
307 i, maxval(alpha), norm_res, t2 - t1
308 ! Convergence in scientific notation
309 !WRITE (iounit, "(T10,I5,T25,1E8.2,T42,1E14.8,T58,F8.2)") &
310 ! i, MAXVAL(alpha), norm_res, t2 - t1
311 CALL m_flush(iounit)
312 END IF
313 END IF
314 IF (converged) THEN
315 IF (iounit > 0 .AND. .NOT. my_silent) THEN
316 WRITE (iounit, "(/,T10,A,I4,A,/)") "The precon solver converged in ", i, " iterations."
317 CALL m_flush(iounit)
318 END IF
319 EXIT iteration
320 END IF
321
322 ! Max number of iteration reached
323 IF (i == max_iter) THEN
324 IF (iounit > 0) THEN
325 WRITE (iounit, "(/,T10,A/)") &
326 "The precon solver didnt converge! Maximum number of iterations reached."
327 CALL m_flush(iounit)
328 END IF
329 converged = .false.
330 END IF
331
332 END DO iteration
333
334 ! Matrix projector
335 CALL projector(qs_env, matrix_p, matrix_cg_z, eps_filter)
336
337 ! Release matrices
338 CALL dbcsr_deallocate_matrix_set(matrix_ax)
339 CALL dbcsr_deallocate_matrix_set(matrix_b)
340 CALL dbcsr_deallocate_matrix_set(matrix_res)
341 CALL dbcsr_deallocate_matrix_set(matrix_cg)
342
343 DEALLOCATE (alpha, beta, new_norm, norm_ca, norm_rr)
344
345 CALL timestop(handle)
346
347 END SUBROUTINE ec_preconditioner
348
349! **************************************************************************************************
350!> \brief AO-based conjugate gradient linear response solver.
351!> In goes the right hand side B of the equation AZ=B, and the linear transformation of the
352!> Hessian matrix A on trial matrices is iteratively solved. Result are
353!> the response density matrix_pz, and the energy-weighted response density matrix_wz.
354!>
355!> \param qs_env ...
356!> \param p_env ...
357!> \param matrix_hz Right hand-side of linear response equation
358!> \param matrix_pz Response density
359!> \param matrix_wz Energy-weighted response density matrix
360!> \param iounit ...
361!> \param should_stop ...
362!>
363!> \param silent ...
364!> \date 01.2020
365!> \author Fabian Belleflamme
366! **************************************************************************************************
367 SUBROUTINE ec_response_ao(qs_env, p_env, matrix_hz, matrix_pz, matrix_wz, iounit, &
368 should_stop, silent)
369
370 TYPE(qs_environment_type), POINTER :: qs_env
371 TYPE(qs_p_env_type), POINTER :: p_env
372 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
373 POINTER :: matrix_hz
374 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
375 POINTER :: matrix_pz, matrix_wz
376 INTEGER, INTENT(IN) :: iounit
377 LOGICAL, INTENT(OUT) :: should_stop
378 LOGICAL, INTENT(IN), OPTIONAL :: silent
379
380 CHARACTER(len=*), PARAMETER :: routinen = 'ec_response_ao'
381
382 INTEGER :: handle, i, ispin, max_iter_lanczos, nao, &
383 nspins, s_sqrt_method, s_sqrt_order
384 LOGICAL :: my_silent, restart
385 REAL(kind=dp) :: eps_filter, eps_lanczos, focc, &
386 min_shift, norm_res, old_conv, shift, &
387 t1, t2
388 REAL(kind=dp), DIMENSION(:), POINTER :: alpha, beta, new_norm, norm_ca, norm_rr, &
389 tr_rz00
390 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_ax, matrix_cg, matrix_cg_z, &
391 matrix_ks, matrix_nsc, matrix_p, matrix_res, matrix_s, matrix_z, matrix_z0, rho_ao
392 TYPE(dbcsr_type) :: matrix_s_sqrt, matrix_s_sqrt_inv, &
393 matrix_tmp
394 TYPE(dft_control_type), POINTER :: dft_control
395 TYPE(linres_control_type), POINTER :: linres_control
396 TYPE(qs_rho_type), POINTER :: rho
397 TYPE(section_vals_type), POINTER :: solver_section
398
399 CALL timeset(routinen, handle)
400
401 my_silent = .false.
402 IF (PRESENT(silent)) my_silent = silent
403
404 cpassert(ASSOCIATED(qs_env))
405 cpassert(ASSOCIATED(matrix_hz))
406 cpassert(ASSOCIATED(matrix_pz))
407 cpassert(ASSOCIATED(matrix_wz))
408
409 NULLIFY (dft_control, ksmat, matrix_s, linres_control, rho)
410
411 t1 = m_walltime()
412
413 CALL get_qs_env(qs_env=qs_env, &
414 dft_control=dft_control, &
415 linres_control=linres_control, &
416 matrix_ks=ksmat, &
417 matrix_s=matrix_s, &
418 rho=rho)
419 nspins = dft_control%nspins
420
421 CALL dbcsr_get_info(matrix_s(1)%matrix, nfullrows_total=nao)
422
423 solver_section => section_vals_get_subs_vals(qs_env%input, "DFT%ENERGY_CORRECTION%RESPONSE_SOLVER")
424 CALL section_vals_val_get(solver_section, "S_SQRT_METHOD", i_val=s_sqrt_method)
425 CALL section_vals_val_get(solver_section, "S_SQRT_ORDER", i_val=s_sqrt_order)
426 CALL section_vals_val_get(solver_section, "EPS_LANCZOS", r_val=eps_lanczos)
427 CALL section_vals_val_get(solver_section, "MAX_ITER_LANCZOS", i_val=max_iter_lanczos)
428
429 eps_filter = linres_control%eps_filter
430
431 CALL qs_rho_get(rho, rho_ao=rho_ao)
432
433 ALLOCATE (alpha(nspins), beta(nspins), new_norm(nspins), norm_ca(nspins), norm_rr(nspins))
434 ALLOCATE (tr_rz00(nspins))
435
436 ! local matrix P, KS, and NSC
437 ! to bring into orthogonal basis
438 NULLIFY (matrix_p, matrix_ks, matrix_nsc)
439 CALL dbcsr_allocate_matrix_set(matrix_p, nspins)
440 CALL dbcsr_allocate_matrix_set(matrix_ks, nspins)
441 CALL dbcsr_allocate_matrix_set(matrix_nsc, nspins)
442 DO ispin = 1, nspins
443 ALLOCATE (matrix_p(ispin)%matrix)
444 ALLOCATE (matrix_ks(ispin)%matrix)
445 ALLOCATE (matrix_nsc(ispin)%matrix)
446 CALL dbcsr_create(matrix_p(ispin)%matrix, name="P_IN ORTHO", &
447 template=ksmat(1)%matrix, &
448 matrix_type=dbcsr_type_no_symmetry)
449 CALL dbcsr_create(matrix_ks(ispin)%matrix, name="KS_IN ORTHO", &
450 template=ksmat(1)%matrix, &
451 matrix_type=dbcsr_type_no_symmetry)
452 CALL dbcsr_create(matrix_nsc(ispin)%matrix, name="NSC IN ORTHO", &
453 template=ksmat(1)%matrix, &
454 matrix_type=dbcsr_type_no_symmetry)
455
456 CALL dbcsr_desymmetrize(rho_ao(ispin)%matrix, matrix_p(ispin)%matrix)
457 CALL dbcsr_desymmetrize(ksmat(ispin)%matrix, matrix_ks(ispin)%matrix)
458 CALL dbcsr_desymmetrize(matrix_hz(ispin)%matrix, matrix_nsc(ispin)%matrix)
459 END DO
460
461 ! Scale matrix_p by factor 1/2 in closed-shell
462 IF (nspins == 1) CALL dbcsr_scale(matrix_p(1)%matrix, 0.5_dp)
463
464 ! Transform P, KS, and Harris kernel matrix into Orthonormal basis
465 CALL dbcsr_create(matrix_s_sqrt, template=matrix_s(1)%matrix, &
466 matrix_type=dbcsr_type_no_symmetry)
467 CALL dbcsr_create(matrix_s_sqrt_inv, template=matrix_s(1)%matrix, &
468 matrix_type=dbcsr_type_no_symmetry)
469
470 SELECT CASE (s_sqrt_method)
471 CASE (ls_s_sqrt_proot)
472 CALL matrix_sqrt_proot(matrix_s_sqrt, matrix_s_sqrt_inv, &
473 matrix_s(1)%matrix, eps_filter, &
474 s_sqrt_order, eps_lanczos, max_iter_lanczos, symmetrize=.true.)
475 CASE (ls_s_sqrt_ns)
476 CALL matrix_sqrt_newton_schulz(matrix_s_sqrt, matrix_s_sqrt_inv, &
477 matrix_s(1)%matrix, eps_filter, &
478 s_sqrt_order, eps_lanczos, max_iter_lanczos)
479 CASE DEFAULT
480 cpabort("Unknown sqrt method.")
481 END SELECT
482
483 ! Transform into orthonormal Lowdin basis
484 DO ispin = 1, nspins
485 CALL transform_m_orth(matrix_p(ispin)%matrix, matrix_s_sqrt, eps_filter)
486 CALL transform_m_orth(matrix_ks(ispin)%matrix, matrix_s_sqrt_inv, eps_filter)
487 CALL transform_m_orth(matrix_nsc(ispin)%matrix, matrix_s_sqrt_inv, eps_filter)
488 END DO
489
490 !----------------------------------------
491 ! Create non-symmetric work matrices: Ax, cg, res
492 ! Content of Ax, cg, cg_z, res, z0 anti-symmetric
493 ! Content of z symmetric
494 !----------------------------------------
495
496 CALL dbcsr_create(matrix_tmp, template=matrix_s(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
497
498 NULLIFY (matrix_ax, matrix_cg, matrix_cg_z, matrix_res, matrix_z, matrix_z0)
499 CALL dbcsr_allocate_matrix_set(matrix_ax, nspins)
500 CALL dbcsr_allocate_matrix_set(matrix_cg, nspins)
501 CALL dbcsr_allocate_matrix_set(matrix_cg_z, nspins)
502 CALL dbcsr_allocate_matrix_set(matrix_res, nspins)
503 CALL dbcsr_allocate_matrix_set(matrix_z, nspins)
504 CALL dbcsr_allocate_matrix_set(matrix_z0, nspins)
505
506 DO ispin = 1, nspins
507 ALLOCATE (matrix_ax(ispin)%matrix)
508 ALLOCATE (matrix_cg(ispin)%matrix)
509 ALLOCATE (matrix_cg_z(ispin)%matrix)
510 ALLOCATE (matrix_res(ispin)%matrix)
511 ALLOCATE (matrix_z(ispin)%matrix)
512 ALLOCATE (matrix_z0(ispin)%matrix)
513 CALL dbcsr_create(matrix_ax(ispin)%matrix, name="linop MATRIX", &
514 template=matrix_s(1)%matrix, &
515 matrix_type=dbcsr_type_no_symmetry)
516 CALL dbcsr_create(matrix_cg(ispin)%matrix, name="TRIAL MATRIX", &
517 template=matrix_s(1)%matrix, &
518 matrix_type=dbcsr_type_no_symmetry)
519 CALL dbcsr_create(matrix_cg_z(ispin)%matrix, name="MATRIX CG-Z", &
520 template=matrix_s(1)%matrix, &
521 matrix_type=dbcsr_type_no_symmetry)
522 CALL dbcsr_create(matrix_res(ispin)%matrix, name="RESIDUE", &
523 template=matrix_s(1)%matrix, &
524 matrix_type=dbcsr_type_no_symmetry)
525 CALL dbcsr_create(matrix_z(ispin)%matrix, name="Z-Matrix", &
526 template=matrix_s(1)%matrix, &
527 matrix_type=dbcsr_type_no_symmetry)
528 CALL dbcsr_create(matrix_z0(ispin)%matrix, name="p after precondi-Matrix", &
529 template=matrix_s(1)%matrix, &
530 matrix_type=dbcsr_type_no_symmetry)
531 END DO
532
533 !----------------------------------------
534 ! Get righ-hand-side operators
535 !----------------------------------------
536
537 ! Spin factor
538 focc = -2.0_dp
539 IF (nspins == 1) focc = -4.0_dp
540
541 ! E^[1]_Harris = -4*G[\delta P]*Pin - Pin*G[\delta P] = -4*[G[\delta P], Pin]
542 CALL commutator(matrix_nsc, matrix_p, matrix_res, eps_filter, .false., alpha=focc)
543
544 ! Initial guess cg_Z
545 DO ispin = 1, nspins
546 CALL dbcsr_copy(matrix_cg_z(ispin)%matrix, matrix_res(ispin)%matrix)
547 END DO
548
549 ! Projector on trial matrix
550 CALL projector(qs_env, matrix_p, matrix_cg_z, eps_filter)
551
552 ! Ax0
553 CALL build_hessian_op(qs_env=qs_env, &
554 p_env=p_env, &
555 matrix_ks=matrix_ks, &
556 matrix_p=matrix_p, & ! p
557 matrix_s_sqrt_inv=matrix_s_sqrt_inv, &
558 matrix_cg=matrix_cg_z, & ! cg
559 matrix_ax=matrix_ax, &
560 eps_filter=eps_filter)
561
562 ! r_0 = b - Ax0
563 DO ispin = 1, nspins
564 CALL dbcsr_add(matrix_res(ispin)%matrix, matrix_ax(ispin)%matrix, 1.0_dp, -1.0_dp)
565 END DO
566
567 ! Matrix projector T
568 CALL projector(qs_env, matrix_p, matrix_res, eps_filter)
569
570 ! Preconditioner
571 linres_control%flag = ""
572 IF (linres_control%preconditioner_type == precond_mlp) THEN
573 ! M * z_0 = r_0
574 ! Conjugate gradient returns z_0
575 CALL ec_preconditioner(qs_env=qs_env, &
576 matrix_ks=matrix_ks, &
577 matrix_p=matrix_p, &
578 matrix_rhs=matrix_res, &
579 matrix_cg_z=matrix_z0, &
580 eps_filter=eps_filter, &
581 iounit=iounit, silent=silent)
582 linres_control%flag = "PCG-AO"
583 ELSE
584 ! z_0 = r_0
585 DO ispin = 1, nspins
586 CALL dbcsr_copy(matrix_z0(ispin)%matrix, matrix_res(ispin)%matrix)
587 linres_control%flag = "CG-AO"
588 END DO
589 END IF
590
591 norm_res = 0.0_dp
592
593 DO ispin = 1, nspins
594 ! cg = p_0 = z_0
595 CALL dbcsr_copy(matrix_cg(ispin)%matrix, matrix_z0(ispin)%matrix)
596
597 ! Tr(r_0 * z_0)
598 CALL dbcsr_dot(matrix_res(ispin)%matrix, matrix_cg(ispin)%matrix, norm_rr(ispin))
599
600 IF (norm_rr(ispin) < 0.0_dp) cpabort("norm_rr < 0")
601 norm_res = max(norm_res, abs(norm_rr(ispin)/real(nao, dp)))
602 END DO
603
604 ! eigenvalue shifting
605 min_shift = 0.0_dp
606 old_conv = norm_rr(1)
607 shift = min(10.0_dp, max(min_shift, 0.05_dp*old_conv))
608 old_conv = 100.0_dp
609
610 ! header
611 IF (iounit > 0 .AND. .NOT. my_silent) THEN
612 WRITE (iounit, "(/,T3,A,T16,A,T25,A,T38,A,T52,A,/,T3,A)") &
613 "Iteration", "Method", "Stepsize", "Convergence", "Time", &
614 repeat("-", 80)
615 END IF
616
617 alpha(:) = 0.0_dp
618 restart = .false.
619 should_stop = .false.
620 linres_control%converged = .false.
621
622 ! start iteration
623 iteration: DO i = 1, linres_control%max_iter
624
625 ! Convergence criteria
626 ! default for eps 10E-6 in MO_linres
627 IF (norm_res < linres_control%eps) THEN
628 linres_control%converged = .true.
629 END IF
630
631 t2 = m_walltime()
632 IF (i == 1 .OR. mod(i, 1) == 0 .OR. linres_control%converged &
633 .OR. restart .OR. should_stop) THEN
634 IF (iounit > 0 .AND. .NOT. my_silent) THEN
635 WRITE (iounit, "(T5,I5,T18,A3,T28,L1,T38,1E8.2,T48,F16.10,T68,F8.2)") &
636 i, linres_control%flag, restart, maxval(alpha), norm_res, t2 - t1
637 CALL m_flush(iounit)
638 END IF
639 END IF
640 IF (linres_control%converged) THEN
641 IF (iounit > 0) THEN
642 WRITE (iounit, "(/,T2,A,I4,A,T73,F8.2,/)") "The linear solver converged in ", &
643 i, " iterations.", t2 - t1
644 CALL m_flush(iounit)
645 END IF
646 EXIT iteration
647 ELSE IF (should_stop) THEN
648 IF (iounit > 0) THEN
649 WRITE (iounit, "(/,T2,A,I4,A,/)") "The linear solver did NOT converge! External stop"
650 CALL m_flush(iounit)
651 END IF
652 EXIT iteration
653 END IF
654
655 ! Max number of iteration reached
656 IF (i == linres_control%max_iter) THEN
657 IF (iounit > 0) THEN
658 WRITE (iounit, "(/,T2,A/)") &
659 "The linear solver didnt converge! Maximum number of iterations reached."
660 CALL m_flush(iounit)
661 END IF
662 linres_control%converged = .false.
663 END IF
664
665 ! Hessian Ax = [F,B] + [G(B),P]
666 CALL build_hessian_op(qs_env=qs_env, &
667 p_env=p_env, &
668 matrix_ks=matrix_ks, &
669 matrix_p=matrix_p, & ! p
670 matrix_s_sqrt_inv=matrix_s_sqrt_inv, &
671 matrix_cg=matrix_cg, & ! cg
672 matrix_ax=matrix_ax, &
673 eps_filter=eps_filter)
674
675 ! Matrix projector T
676 CALL projector(qs_env, matrix_p, matrix_ax, eps_filter)
677
678 DO ispin = 1, nspins
679
680 CALL dbcsr_filter(matrix_ax(ispin)%matrix, eps_filter)
681 ! norm_cA = tr(Ap_j * p_j)
682 CALL dbcsr_dot(matrix_cg(ispin)%matrix, matrix_ax(ispin)%matrix, norm_ca(ispin))
683
684 IF (norm_ca(ispin) < 0.0_dp) THEN
685
686 ! Recalculate w/o preconditioner
687 IF (i > 1) THEN
688 ! p = -z + beta*p
689 CALL dbcsr_add(matrix_cg(ispin)%matrix, matrix_z0(ispin)%matrix, &
690 beta(ispin), -1.0_dp)
691 CALL dbcsr_dot(matrix_res(ispin)%matrix, matrix_res(ispin)%matrix, new_norm(ispin))
692 beta(ispin) = new_norm(ispin)/tr_rz00(ispin)
693 CALL dbcsr_add(matrix_cg(ispin)%matrix, matrix_res(ispin)%matrix, &
694 beta(ispin), 1.0_dp)
695 norm_rr(ispin) = new_norm(ispin)
696 ELSE
697 CALL dbcsr_copy(matrix_res(ispin)%matrix, matrix_cg(ispin)%matrix)
698 CALL dbcsr_dot(matrix_res(ispin)%matrix, matrix_res(ispin)%matrix, norm_rr(ispin))
699 END IF
700
701 CALL build_hessian_op(qs_env=qs_env, &
702 p_env=p_env, &
703 matrix_ks=matrix_ks, &
704 matrix_p=matrix_p, & ! p
705 matrix_s_sqrt_inv=matrix_s_sqrt_inv, &
706 matrix_cg=matrix_cg, & ! cg
707 matrix_ax=matrix_ax, &
708 eps_filter=eps_filter)
709
710 ! Matrix projector T
711 CALL projector(qs_env, matrix_p, matrix_ax, eps_filter)
712
713 CALL dbcsr_dot(matrix_cg(ispin)%matrix, matrix_ax(ispin)%matrix, norm_ca(ispin))
714
715 cpabort("tr(Ap_j*p_j) < 0")
716 IF (abnormal_value(norm_ca(ispin))) THEN
717 cpabort("Preconditioner: Tr[Ap_j*p_j] is an abnormal value (NaN/Inf)")
718 END IF
719
720 END IF
721
722 END DO
723
724 DO ispin = 1, nspins
725 ! Determine step-size
726 IF (norm_ca(ispin) < linres_control%eps) THEN
727 alpha(ispin) = 1.0_dp
728 ELSE
729 alpha(ispin) = norm_rr(ispin)/norm_ca(ispin)
730 END IF
731
732 ! x_j+1 = x_j + alpha*p_j
733 ! save response-denisty of this iteration
734 CALL dbcsr_add(matrix_cg_z(ispin)%matrix, matrix_cg(ispin)%matrix, 1.0_dp, alpha(ispin))
735 END DO
736
737 ! need to recompute the residue
738 restart = .false.
739 IF (mod(i, linres_control%restart_every) == 0) THEN
740 !
741 ! r_j+1 = b - A * x_j+1
742 CALL build_hessian_op(qs_env=qs_env, &
743 p_env=p_env, &
744 matrix_ks=matrix_ks, &
745 matrix_p=matrix_p, &
746 matrix_s_sqrt_inv=matrix_s_sqrt_inv, &
747 matrix_cg=matrix_cg_z, & ! cg
748 matrix_ax=matrix_ax, &
749 eps_filter=eps_filter)
750 ! b
751 CALL commutator(matrix_nsc, matrix_p, matrix_res, eps_filter, .false., alpha=focc)
752
753 DO ispin = 1, nspins
754 CALL dbcsr_add(matrix_res(ispin)%matrix, matrix_ax(ispin)%matrix, 1.0_dp, -1.0_dp)
755 END DO
756
757 CALL projector(qs_env, matrix_p, matrix_res, eps_filter)
758 !
759 restart = .true.
760 ELSE
761 ! proj Ap onto the virtual subspace
762 CALL projector(qs_env, matrix_p, matrix_ax, eps_filter)
763 !
764 ! r_j+1 = r_j - alpha * Ap_j
765 DO ispin = 1, nspins
766 CALL dbcsr_add(matrix_res(ispin)%matrix, matrix_ax(ispin)%matrix, 1.0_dp, -alpha(ispin))
767 END DO
768 restart = .false.
769 END IF
770
771 ! Preconditioner
772 linres_control%flag = ""
773 IF (linres_control%preconditioner_type == precond_mlp) THEN
774 ! M * z_j+1 = r_j+1
775 ! Conjugate gradient returns z_j+1
776 CALL ec_preconditioner(qs_env=qs_env, &
777 matrix_ks=matrix_ks, &
778 matrix_p=matrix_p, &
779 matrix_rhs=matrix_res, &
780 matrix_cg_z=matrix_z0, &
781 eps_filter=eps_filter, &
782 iounit=iounit, silent=silent)
783 linres_control%flag = "PCG-AO"
784 ELSE
785 DO ispin = 1, nspins
786 CALL dbcsr_copy(matrix_z0(ispin)%matrix, matrix_res(ispin)%matrix)
787 END DO
788 linres_control%flag = "CG-AO"
789 END IF
790
791 norm_res = 0.0_dp
792
793 DO ispin = 1, nspins
794 ! Tr[r_j+1*z_j+1]
795 CALL dbcsr_dot(matrix_res(ispin)%matrix, matrix_z0(ispin)%matrix, new_norm(ispin))
796 IF (new_norm(ispin) < 0.0_dp) cpabort("tr(r_j+1*z_j+1) < 0")
797 IF (abnormal_value(new_norm(ispin))) THEN
798 cpabort("Preconditioner: Tr[r_j+1*z_j+1] is an abnormal value (NaN/Inf)")
799 END IF
800 norm_res = max(norm_res, new_norm(ispin)/real(nao, dp))
801
802 IF (norm_rr(ispin) < linres_control%eps .OR. new_norm(ispin) < linres_control%eps) THEN
803 beta(ispin) = 0.0_dp
804 linres_control%converged = .true.
805 ELSE
806 beta(ispin) = new_norm(ispin)/norm_rr(ispin)
807 END IF
808
809 ! update new search vector (matrix cg)
810 ! Here: cg_j+1 = z_j+1 + beta*cg_j
811 CALL dbcsr_add(matrix_cg(ispin)%matrix, matrix_z0(ispin)%matrix, beta(ispin), 1.0_dp)
812 CALL dbcsr_filter(matrix_cg(ispin)%matrix, eps_filter)
813
814 tr_rz00(ispin) = norm_rr(ispin)
815 norm_rr(ispin) = new_norm(ispin)
816 END DO
817
818 ! Can we exit the loop?
819 CALL external_control(should_stop, "LS_SOLVER", target_time=qs_env%target_time, &
820 start_time=qs_env%start_time)
821
822 END DO iteration
823
824 ! Matrix projector
825 CALL projector(qs_env, matrix_p, matrix_cg_z, eps_filter)
826
827 ! Z = [cg_z,P]
828 CALL commutator(matrix_cg_z, matrix_p, matrix_z, eps_filter, .true., alpha=0.5_dp)
829
830 DO ispin = 1, nspins
831 ! Transform Z-matrix back into non-orthogonal basis
832 CALL transform_m_orth(matrix_z(ispin)%matrix, matrix_s_sqrt_inv, eps_filter)
833
834 ! Export Z-Matrix
835 CALL dbcsr_copy(matrix_pz(ispin)%matrix, matrix_z(ispin)%matrix, keep_sparsity=.true.)
836 END DO
837
838 ! Calculate energy-weighted response density matrix
839 ! AO: Wz = 0.5*(Z*KS*P + P*KS*Z)
840 CALL ec_wz_matrix(qs_env, matrix_pz, matrix_wz, eps_filter)
841
842 ! Release matrices
843 CALL dbcsr_release(matrix_tmp)
844
845 CALL dbcsr_release(matrix_s_sqrt)
846 CALL dbcsr_release(matrix_s_sqrt_inv)
847
848 CALL dbcsr_deallocate_matrix_set(matrix_p)
849 CALL dbcsr_deallocate_matrix_set(matrix_ks)
850 CALL dbcsr_deallocate_matrix_set(matrix_nsc)
851 CALL dbcsr_deallocate_matrix_set(matrix_z)
852 CALL dbcsr_deallocate_matrix_set(matrix_ax)
853 CALL dbcsr_deallocate_matrix_set(matrix_res)
854 CALL dbcsr_deallocate_matrix_set(matrix_cg)
855 CALL dbcsr_deallocate_matrix_set(matrix_cg_z)
856 CALL dbcsr_deallocate_matrix_set(matrix_z0)
857
858 DEALLOCATE (alpha, beta, new_norm, norm_ca, norm_rr)
859 DEALLOCATE (tr_rz00)
860
861 CALL timestop(handle)
862
863 END SUBROUTINE ec_response_ao
864
865! **************************************************************************************************
866!> \brief Compute matrix_wz as needed for the forces
867!> Wz = 0.5*(Z*KS*P + P*KS*Z) (closed-shell)
868!> \param qs_env ...
869!> \param matrix_z The response density we just calculated
870!> \param matrix_wz The energy weighted response-density matrix
871!> \param eps_filter ...
872!> \par History
873!> 2020.2 created [Fabian Belleflamme]
874!> \author Fabian Belleflamme
875! **************************************************************************************************
876 SUBROUTINE ec_wz_matrix(qs_env, matrix_z, matrix_wz, eps_filter)
877
878 TYPE(qs_environment_type), POINTER :: qs_env
879 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
880 POINTER :: matrix_z
881 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
882 POINTER :: matrix_wz
883 REAL(kind=dp), INTENT(IN) :: eps_filter
884
885 CHARACTER(len=*), PARAMETER :: routinen = 'ec_wz_matrix'
886
887 INTEGER :: handle, ispin, nspins
888 REAL(kind=dp) :: scaling
889 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p, matrix_s
890 TYPE(dbcsr_type) :: matrix_tmp, matrix_tmp2
891 TYPE(dft_control_type), POINTER :: dft_control
892 TYPE(qs_rho_type), POINTER :: rho
893
894 CALL timeset(routinen, handle)
895
896 cpassert(ASSOCIATED(qs_env))
897 cpassert(ASSOCIATED(matrix_z))
898 cpassert(ASSOCIATED(matrix_wz))
899
900 CALL get_qs_env(qs_env=qs_env, &
901 dft_control=dft_control, &
902 matrix_ks=matrix_ks, &
903 matrix_s=matrix_s, &
904 rho=rho)
905 nspins = dft_control%nspins
906
907 CALL qs_rho_get(rho, rho_ao=matrix_p)
908
909 ! Init temp matrices
910 CALL dbcsr_create(matrix_tmp, template=matrix_z(1)%matrix, &
911 matrix_type=dbcsr_type_no_symmetry)
912 CALL dbcsr_create(matrix_tmp2, template=matrix_z(1)%matrix, &
913 matrix_type=dbcsr_type_no_symmetry)
914
915 ! Scale matrix_p by factor 1/2 in closed-shell
916 scaling = 1.0_dp
917 IF (nspins == 1) scaling = 0.5_dp
918
919 ! Whz = ZFP + PFZ = Z(FP) + (Z(FP))^T
920 DO ispin = 1, nspins
921
922 ! tmp = FP
923 CALL dbcsr_multiply("N", "N", scaling, matrix_ks(ispin)%matrix, matrix_p(ispin)%matrix, &
924 0.0_dp, matrix_tmp, filter_eps=eps_filter, retain_sparsity=.false.)
925
926 ! tmp2 = ZFP
927 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_z(ispin)%matrix, matrix_tmp, &
928 0.0_dp, matrix_tmp2, filter_eps=eps_filter, retain_sparsity=.false.)
929
930 ! tmp = (ZFP)^T
931 CALL dbcsr_transposed(matrix_tmp, matrix_tmp2)
932
933 ! tmp = ZFP + (ZFP)^T
934 CALL dbcsr_add(matrix_tmp, matrix_tmp2, 1.0_dp, 1.0_dp)
935
936 CALL dbcsr_filter(matrix_tmp, eps_filter)
937
938 ! Whz = ZFP + PFZ
939 CALL dbcsr_copy(matrix_wz(ispin)%matrix, matrix_tmp, keep_sparsity=.true.)
940
941 END DO
942
943 ! Release matrices
944 CALL dbcsr_release(matrix_tmp)
945 CALL dbcsr_release(matrix_tmp2)
946
947 CALL timestop(handle)
948
949 END SUBROUTINE ec_wz_matrix
950
951! **************************************************************************************************
952!> \brief Calculate first term of electronic Hessian M = [F, B]
953!> acting as liner transformation on trial matrix (matrix_cg)
954!> with intermediate response density B = [cg,P] = cg*P - P*cg = cg*P + (cg*P)^T
955!>
956!> All matrices are in orthonormal basis
957!>
958!> \param matrix_ks Ground-state Kohn-Sham matrix
959!> \param matrix_p Ground-state Density matrix
960!> \param matrix_cg Trial matrix
961!> \param matrix_b Intermediate response density
962!> \param matrix_Ax First term of electronic Hessian applied on trial matrix (matrix_cg)
963!>
964!> \param eps_filter ...
965!> \date 12.2019
966!> \author Fabian Belleflamme
967! **************************************************************************************************
968 SUBROUTINE hessian_op1(matrix_ks, matrix_p, matrix_cg, matrix_b, matrix_Ax, eps_filter)
969
970 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
971 POINTER :: matrix_ks, matrix_p, matrix_cg
972 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
973 POINTER :: matrix_b, matrix_ax
974 REAL(kind=dp), INTENT(IN) :: eps_filter
975
976 CHARACTER(len=*), PARAMETER :: routinen = 'hessian_op1'
977
978 INTEGER :: handle
979
980 CALL timeset(routinen, handle)
981
982 cpassert(ASSOCIATED(matrix_ks))
983 cpassert(ASSOCIATED(matrix_p))
984 cpassert(ASSOCIATED(matrix_cg))
985 cpassert(ASSOCIATED(matrix_b))
986 cpassert(ASSOCIATED(matrix_ax))
987
988 ! Build intermediate density matrix
989 ! B = [cg, P] = cg*P - P*cg = cg*P + (cg*P)^T
990 CALL commutator(matrix_cg, matrix_p, matrix_b, eps_filter, .true.)
991
992 ! Build first part of operator
993 ! Ax = [F,[cg,P]] = [F,B]
994 CALL commutator(matrix_ks, matrix_b, matrix_ax, eps_filter, .false.)
995
996 CALL timestop(handle)
997
998 END SUBROUTINE hessian_op1
999
1000! **************************************************************************************************
1001!> \brief calculate linear transformation of Hessian matrix on a trial matrix matrix_cg
1002!> which is stored in response density B = [cg,P] = cg*P - P*cg = cg*P + (cg*P)^T
1003!> Ax = [F, B] + [G(B), Pin] in orthonormal basis
1004!>
1005!> \param qs_env ...
1006!> \param p_env ...
1007!> \param matrix_ks Ground-state Kohn-Sham matrix
1008!> \param matrix_p Ground-state Density matrix
1009!> \param matrix_s_sqrt_inv S^(-1/2) needed for transformation to/from orthonormal basis
1010!> \param matrix_cg Trial matrix
1011!> \param matrix_Ax Electronic Hessian applied on trial matrix (matrix_cg)
1012!> \param eps_filter ...
1013!>
1014!> \date 12.2019
1015!> \author Fabian Belleflamme
1016! **************************************************************************************************
1017 SUBROUTINE build_hessian_op(qs_env, p_env, matrix_ks, matrix_p, matrix_s_sqrt_inv, &
1018 matrix_cg, matrix_Ax, eps_filter)
1019
1020 TYPE(qs_environment_type), POINTER :: qs_env
1021 TYPE(qs_p_env_type), POINTER :: p_env
1022 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1023 POINTER :: matrix_ks, matrix_p
1024 TYPE(dbcsr_type), INTENT(IN) :: matrix_s_sqrt_inv
1025 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1026 POINTER :: matrix_cg
1027 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
1028 POINTER :: matrix_ax
1029 REAL(kind=dp), INTENT(IN) :: eps_filter
1030
1031 CHARACTER(len=*), PARAMETER :: routinen = 'build_hessian_op'
1032
1033 INTEGER :: handle, ispin, nspins
1034 REAL(kind=dp) :: chksum
1035 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_b, rho1_ao
1036 TYPE(dft_control_type), POINTER :: dft_control
1037 TYPE(mp_para_env_type), POINTER :: para_env
1038 TYPE(qs_rho_type), POINTER :: rho
1039
1040 CALL timeset(routinen, handle)
1041
1042 cpassert(ASSOCIATED(qs_env))
1043 cpassert(ASSOCIATED(matrix_ks))
1044 cpassert(ASSOCIATED(matrix_p))
1045 cpassert(ASSOCIATED(matrix_cg))
1046 cpassert(ASSOCIATED(matrix_ax))
1047
1048 CALL get_qs_env(qs_env=qs_env, &
1049 dft_control=dft_control, &
1050 para_env=para_env, &
1051 rho=rho)
1052 nspins = dft_control%nspins
1053
1054 NULLIFY (matrix_b)
1055 CALL dbcsr_allocate_matrix_set(matrix_b, nspins)
1056 DO ispin = 1, nspins
1057 ALLOCATE (matrix_b(ispin)%matrix)
1058 CALL dbcsr_create(matrix_b(ispin)%matrix, name="[X,P] RSP DNSTY", &
1059 template=matrix_p(1)%matrix, &
1060 matrix_type=dbcsr_type_no_symmetry)
1061 END DO
1062
1063 ! Build uncoupled term of Hessian linear transformation
1064 CALL hessian_op1(matrix_ks, matrix_p, matrix_cg, matrix_b, matrix_ax, eps_filter)
1065
1066 ! Avoid the buildup of noisy blocks
1067 DO ispin = 1, nspins
1068 CALL dbcsr_filter(matrix_b(ispin)%matrix, eps_filter)
1069 END DO
1070
1071 chksum = 0.0_dp
1072 DO ispin = 1, nspins
1073 chksum = chksum + dbcsr_checksum(matrix_b(ispin)%matrix)
1074 END DO
1075
1076 ! skip the kernel if the DM is very small
1077 IF (chksum > 1.0e-14_dp) THEN
1078
1079 ! Bring matrix B as density on grid
1080
1081 ! prepare perturbation environment
1082 CALL p_env_check_i_alloc(p_env, qs_env)
1083
1084 ! Get response density matrix
1085 CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
1086
1087 DO ispin = 1, nspins
1088 ! Transform B into NON-ortho basis for collocation
1089 CALL transform_m_orth(matrix_b(ispin)%matrix, matrix_s_sqrt_inv, eps_filter)
1090 ! Filter
1091 CALL dbcsr_filter(matrix_b(ispin)%matrix, eps_filter)
1092 ! Keep symmetry of density matrix
1093 CALL dbcsr_copy(rho1_ao(ispin)%matrix, matrix_b(ispin)%matrix, keep_sparsity=.true.)
1094 CALL dbcsr_copy(p_env%p1(ispin)%matrix, matrix_b(ispin)%matrix, keep_sparsity=.true.)
1095 END DO
1096
1097 ! Updates densities on grid wrt density matrix
1098 CALL p_env_update_rho(p_env, qs_env)
1099
1100 DO ispin = 1, nspins
1101 CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
1102 IF (ASSOCIATED(p_env%kpp1_admm)) CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
1103 END DO
1104
1105 ! Calculate kernel
1106 ! Ax = F*B - B*F + G(B)*P - P*G(B)
1107 ! IN/OUT IN IN IN
1108 CALL hessian_op2(qs_env, p_env, matrix_ax, matrix_p, matrix_s_sqrt_inv, eps_filter)
1109
1110 END IF
1111
1112 CALL dbcsr_deallocate_matrix_set(matrix_b)
1113
1114 CALL timestop(handle)
1115
1116 END SUBROUTINE build_hessian_op
1117
1118! **************************************************************************************************
1119!> \brief Calculate lin transformation of Hessian matrix on a trial matrix matrix_cg
1120!> which is stored in response density B = [cg,P] = cg*P - P*cg = cg*P + (cg*P)^T
1121!> Ax = [F, B] + [G(B), Pin] in orthonormal basis
1122!>
1123!> \param qs_env ...
1124!> \param p_env p-environment with trial density environment
1125!> \param matrix_Ax contains first part of Hessian linear transformation, kernel contribution
1126!> is calculated and added in this routine
1127!> \param matrix_p Density matrix in orthogonal basis
1128!> \param matrix_s_sqrt_inv contains matrix S^(-1/2) for switching to orthonormal Lowdin basis
1129!> \param eps_filter ...
1130!>
1131!> \date 12.2019
1132!> \author Fabian Belleflamme
1133! **************************************************************************************************
1134 SUBROUTINE hessian_op2(qs_env, p_env, matrix_Ax, matrix_p, matrix_s_sqrt_inv, eps_filter)
1135
1136 TYPE(qs_environment_type), POINTER :: qs_env
1137 TYPE(qs_p_env_type), POINTER :: p_env
1138 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
1139 POINTER :: matrix_ax
1140 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1141 POINTER :: matrix_p
1142 TYPE(dbcsr_type), INTENT(IN) :: matrix_s_sqrt_inv
1143 REAL(kind=dp), INTENT(IN) :: eps_filter
1144
1145 CHARACTER(len=*), PARAMETER :: routinen = 'hessian_op2'
1146
1147 INTEGER :: handle, ispin, nspins
1148 REAL(kind=dp) :: ekin_mol
1149 TYPE(admm_type), POINTER :: admm_env
1150 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_g, matrix_s, rho1_ao, rho_ao
1151 TYPE(dft_control_type), POINTER :: dft_control
1152 TYPE(mp_para_env_type), POINTER :: para_env
1153 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
1154 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g
1155 TYPE(pw_env_type), POINTER :: pw_env
1156 TYPE(pw_poisson_type), POINTER :: poisson_env
1157 TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
1158 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1159 TYPE(pw_r3d_rs_type) :: v_hartree_rspace
1160 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho_r, tau1_r, v_xc, v_xc_tau
1161 TYPE(pw_r3d_rs_type), POINTER :: weights
1162 TYPE(qs_kpp1_env_type), POINTER :: kpp1_env
1163 TYPE(qs_rho_type), POINTER :: rho, rho_aux
1164 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set
1165 TYPE(section_vals_type), POINTER :: input, xc_section, xc_section_aux
1166
1167 CALL timeset(routinen, handle)
1168
1169 NULLIFY (admm_env, dft_control, input, matrix_s, para_env, rho, rho_r, rho1_g, rho1_r)
1170
1171 CALL get_qs_env(qs_env=qs_env, &
1172 admm_env=admm_env, &
1173 dft_control=dft_control, &
1174 input=input, &
1175 matrix_s=matrix_s, &
1176 para_env=para_env, &
1177 rho=rho)
1178 nspins = dft_control%nspins
1179
1180 cpassert(ASSOCIATED(p_env%kpp1))
1181 cpassert(ASSOCIATED(p_env%kpp1_env))
1182 kpp1_env => p_env%kpp1_env
1183
1184 ! Get non-ortho input density matrix on grid
1185 CALL qs_rho_get(rho, rho_ao=rho_ao)
1186 ! Get non-ortho trial density stored in p_env
1187 CALL qs_rho_get(p_env%rho1, rho_g=rho1_g, rho_r=rho1_r, tau_r=tau1_r)
1188
1189 NULLIFY (pw_env)
1190 CALL get_qs_env(qs_env, pw_env=pw_env)
1191 cpassert(ASSOCIATED(pw_env))
1192
1193 NULLIFY (weights)
1194 CALL get_qs_env(qs_env, xcint_weights=weights)
1195
1196 NULLIFY (auxbas_pw_pool, poisson_env, pw_pools)
1197 ! gets the tmp grids
1198 CALL pw_env_get(pw_env=pw_env, &
1199 auxbas_pw_pool=auxbas_pw_pool, &
1200 pw_pools=pw_pools, &
1201 poisson_env=poisson_env)
1202
1203 ! Calculate the NSC Hartree potential
1204 CALL auxbas_pw_pool%create_pw(pw=v_hartree_gspace)
1205 CALL auxbas_pw_pool%create_pw(pw=rho_tot_gspace)
1206 CALL auxbas_pw_pool%create_pw(pw=v_hartree_rspace)
1207
1208 ! XC-Kernel
1209 NULLIFY (v_xc, v_xc_tau, xc_section)
1210
1211 IF (dft_control%do_admm) THEN
1212 xc_section => admm_env%xc_section_primary
1213 ELSE
1214 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1215 END IF
1216
1217 ! add xc-kernel
1218 CALL qs_fxc_create(qs_env, rho, p_env%rho1, rho0_atom_set, xc_section, .false., &
1219 v_xc, v_xc_tau, rho1_atom_set)
1220
1221 DO ispin = 1, nspins
1222 CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
1223 IF (ASSOCIATED(v_xc_tau)) THEN
1224 CALL pw_scale(v_xc_tau(ispin), v_xc_tau(ispin)%pw_grid%dvol)
1225 END IF
1226 END DO
1227
1228 ! ADMM Correction
1229 IF (dft_control%do_admm) THEN
1230 IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
1231 IF (.NOT. ASSOCIATED(kpp1_env%deriv_set_admm)) THEN
1232 xc_section_aux => admm_env%xc_section_aux
1233 CALL get_admm_env(admm_env, rho_aux_fit=rho_aux)
1234 ALLOCATE (kpp1_env%deriv_set_admm, kpp1_env%rho_set_admm)
1235 CALL qs_fxc_prep(qs_env, rho_aux, kpp1_env%rho_set_admm, kpp1_env%deriv_set_admm, &
1236 xc_section_aux, pw_env, is_triplet=.false.)
1237 END IF
1238 END IF
1239 END IF
1240
1241 ! take trial density to build G^{H}[B]
1242 CALL pw_zero(rho_tot_gspace)
1243 DO ispin = 1, nspins
1244 CALL pw_axpy(rho1_g(ispin), rho_tot_gspace)
1245 END DO
1246
1247 ! get Hartree potential from rho_tot_gspace
1248 CALL pw_poisson_solve(poisson_env, rho_tot_gspace, &
1249 vhartree=v_hartree_gspace)
1250 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
1251 CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
1252
1253 ! Add v_xc + v_H
1254 DO ispin = 1, nspins
1255 CALL pw_axpy(v_hartree_rspace, v_xc(ispin))
1256 END DO
1257 IF (nspins == 1) THEN
1258 CALL pw_scale(v_xc(1), 2.0_dp)
1259 IF (ASSOCIATED(v_xc_tau)) CALL pw_scale(v_xc_tau(1), 2.0_dp)
1260 END IF
1261
1262 DO ispin = 1, nspins
1263 ! Integrate with ground-state density matrix, in non-orthogonal basis
1264 CALL integrate_v_rspace(v_rspace=v_xc(ispin), &
1265 pmat=rho_ao(ispin), &
1266 hmat=p_env%kpp1(ispin), &
1267 qs_env=qs_env, &
1268 calculate_forces=.false., &
1269 basis_type="ORB")
1270 IF (ASSOCIATED(v_xc_tau)) THEN
1271 CALL integrate_v_rspace(v_rspace=v_xc_tau(ispin), &
1272 pmat=rho_ao(ispin), &
1273 hmat=p_env%kpp1(ispin), &
1274 qs_env=qs_env, &
1275 compute_tau=.true., &
1276 calculate_forces=.false., &
1277 basis_type="ORB")
1278 END IF
1279 END DO
1280
1281 ! Hartree-Fock contribution
1282 CALL apply_hfx(qs_env, p_env)
1283 ! Calculate ADMM exchange correction to kernel
1284 CALL apply_xc_admm(qs_env, p_env)
1285 ! Add contribution from ADMM exchange correction to kernel
1286 CALL p_env_finish_kpp1(qs_env, p_env)
1287
1288 ! Calculate KG correction to kernel
1289 IF (dft_control%qs_control%do_kg) THEN
1290 IF (qs_env%kg_env%tnadd_method == kg_tnadd_embed .OR. &
1291 qs_env%kg_env%tnadd_method == kg_tnadd_embed_ri) THEN
1292
1293 cpassert(dft_control%nimages == 1)
1294 ekin_mol = 0.0_dp
1295 CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
1296 CALL kg_ekin_subset(qs_env=qs_env, &
1297 ks_matrix=p_env%kpp1, &
1298 ekin_mol=ekin_mol, &
1299 calc_force=.false., &
1300 do_kernel=.true., &
1301 pmat_ext=rho1_ao)
1302 END IF
1303 END IF
1304
1305 ! Init response kernel matrix
1306 ! matrix G(B)
1307 NULLIFY (matrix_g)
1308 CALL dbcsr_allocate_matrix_set(matrix_g, nspins)
1309 DO ispin = 1, nspins
1310 ALLOCATE (matrix_g(ispin)%matrix)
1311 CALL dbcsr_copy(matrix_g(ispin)%matrix, p_env%kpp1(ispin)%matrix, &
1312 name="MATRIX Kernel")
1313 END DO
1314
1315 ! Transforming G(B) into orthonormal basis
1316 ! Careful, this de-symmetrizes matrix_G
1317 DO ispin = 1, nspins
1318 CALL transform_m_orth(matrix_g(ispin)%matrix, matrix_s_sqrt_inv, eps_filter)
1319 CALL dbcsr_filter(matrix_g(ispin)%matrix, eps_filter)
1320 END DO
1321
1322 ! Hessian already contains Ax = [F,B] (ORTHO), now adding
1323 ! Ax = Ax + G(B)P - (G(B)P)^T
1324 CALL commutator(matrix_g, matrix_p, matrix_ax, eps_filter, .false., 1.0_dp, 1.0_dp)
1325
1326 ! release pw grids
1327 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
1328 CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
1329 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
1330 DO ispin = 1, nspins
1331 CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
1332 END DO
1333 DEALLOCATE (v_xc)
1334 IF (ASSOCIATED(v_xc_tau)) THEN
1335 DO ispin = 1, nspins
1336 CALL auxbas_pw_pool%give_back_pw(v_xc_tau(ispin))
1337 END DO
1338 DEALLOCATE (v_xc_tau)
1339 END IF
1340
1341 CALL dbcsr_deallocate_matrix_set(matrix_g)
1342
1343 CALL timestop(handle)
1344
1345 END SUBROUTINE hessian_op2
1346
1347! **************************************************************************************************
1348!> \brief computes (anti-)commutator exploiting (anti-)symmetry:
1349!> A symmetric : RES = beta*RES + k*[A,B] = k*(AB-(AB)^T)
1350!> A anti-sym : RES = beta*RES + k*{A,B} = k*(AB+(AB)^T)
1351!>
1352!> \param a Matrix A
1353!> \param b Matrix B
1354!> \param res Commutator result
1355!> \param eps_filter filtering threshold for sparse matrices
1356!> \param anticomm Calculate anticommutator
1357!> \param alpha Scaling of anti-/commutator
1358!> \param beta Scaling of inital content of result matrix
1359!>
1360!> \par History
1361!> 2020.07 Fabian Belleflamme (based on commutator_symm)
1362! **************************************************************************************************
1363 SUBROUTINE commutator(a, b, res, eps_filter, anticomm, alpha, beta)
1364
1365 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1366 POINTER :: a, b, res
1367 REAL(kind=dp) :: eps_filter
1368 LOGICAL :: anticomm
1369 REAL(kind=dp), OPTIONAL :: alpha, beta
1370
1371 CHARACTER(LEN=*), PARAMETER :: routinen = 'commutator'
1372
1373 INTEGER :: handle, ispin
1374 REAL(kind=dp) :: facc, myalpha, mybeta
1375 TYPE(dbcsr_type) :: work, work2
1376
1377 CALL timeset(routinen, handle)
1378
1379 cpassert(ASSOCIATED(a))
1380 cpassert(ASSOCIATED(b))
1381 cpassert(ASSOCIATED(res))
1382
1383 CALL dbcsr_create(work, template=a(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1384 CALL dbcsr_create(work2, template=a(1)%matrix, matrix_type=dbcsr_type_no_symmetry)
1385
1386 ! Scaling of anti-/commutator
1387 myalpha = 1.0_dp
1388 IF (PRESENT(alpha)) myalpha = alpha
1389 ! Scaling of result matrix
1390 mybeta = 0.0_dp
1391 IF (PRESENT(beta)) mybeta = beta
1392 ! Add/subtract second term when calculating anti-/commutator
1393 facc = -1.0_dp
1394 IF (anticomm) facc = 1.0_dp
1395
1396 DO ispin = 1, SIZE(a)
1397
1398 CALL dbcsr_multiply("N", "N", myalpha, a(ispin)%matrix, b(ispin)%matrix, &
1399 0.0_dp, work, filter_eps=eps_filter)
1400 CALL dbcsr_transposed(work2, work)
1401
1402 ! RES= beta*RES + alpha*{A,B} = beta*RES + alpha*[AB+(AB)T]
1403 ! RES= beta*RES + alpha*[A,B] = beta*RES + alpha*[AB-(AB)T]
1404 CALL dbcsr_add(work, work2, 1.0_dp, facc)
1405
1406 CALL dbcsr_add(res(ispin)%matrix, work, mybeta, 1.0_dp)
1407
1408 END DO
1409
1410 CALL dbcsr_release(work)
1411 CALL dbcsr_release(work2)
1412
1413 CALL timestop(handle)
1414
1415 END SUBROUTINE commutator
1416
1417! **************************************************************************************************
1418!> \brief Projector P(M) = P*M*Q^T + Q*M*P^T
1419!> with P = D
1420!> with Q = (1-D)
1421!>
1422!> \param qs_env ...
1423!> \param matrix_p Ground-state density in orthonormal basis
1424!> \param matrix_io Matrix to which projector is applied.
1425!>
1426!> \param eps_filter ...
1427!> \date 06.2020
1428!> \author Fabian Belleflamme
1429! **************************************************************************************************
1430 SUBROUTINE projector(qs_env, matrix_p, matrix_io, eps_filter)
1431
1432 TYPE(qs_environment_type), POINTER :: qs_env
1433 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(IN), &
1434 POINTER :: matrix_p
1435 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
1436 POINTER :: matrix_io
1437 REAL(kind=dp), INTENT(IN) :: eps_filter
1438
1439 CHARACTER(len=*), PARAMETER :: routinen = 'projector'
1440
1441 INTEGER :: handle, ispin, nspins
1442 TYPE(dbcsr_type) :: matrix_q, matrix_tmp
1443 TYPE(dft_control_type), POINTER :: dft_control
1444 TYPE(mp_para_env_type), POINTER :: para_env
1445
1446 CALL timeset(routinen, handle)
1447
1448 CALL get_qs_env(qs_env=qs_env, &
1449 dft_control=dft_control, &
1450 para_env=para_env)
1451 nspins = dft_control%nspins
1452
1453 CALL dbcsr_create(matrix_q, template=matrix_p(1)%matrix, &
1454 matrix_type=dbcsr_type_no_symmetry)
1455 CALL dbcsr_create(matrix_tmp, template=matrix_p(1)%matrix, &
1456 matrix_type=dbcsr_type_no_symmetry)
1457
1458 ! Q = (1 - P)
1459 CALL dbcsr_copy(matrix_q, matrix_p(1)%matrix)
1460 CALL dbcsr_scale(matrix_q, -1.0_dp)
1461 CALL dbcsr_add_on_diag(matrix_q, 1.0_dp)
1462 CALL dbcsr_finalize(matrix_q)
1463
1464 ! Proj(M) = P*M*Q + Q*M*P
1465 ! with P = D = CC^T
1466 ! and Q = (1 - P)
1467 DO ispin = 1, nspins
1468
1469 ! tmp1 = P*M
1470 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_p(ispin)%matrix, matrix_io(ispin)%matrix, &
1471 0.0_dp, matrix_tmp, filter_eps=eps_filter)
1472 ! m_io = P*M*Q
1473 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_tmp, matrix_q, &
1474 0.0_dp, matrix_io(ispin)%matrix, filter_eps=eps_filter)
1475
1476 ! tmp = (P^T*M^T*Q^T)^T = -(P*M*Q)^T
1477 CALL dbcsr_transposed(matrix_tmp, matrix_io(ispin)%matrix)
1478 CALL dbcsr_add(matrix_io(ispin)%matrix, matrix_tmp, 1.0_dp, -1.0_dp)
1479
1480 END DO
1481
1482 CALL dbcsr_release(matrix_tmp)
1483 CALL dbcsr_release(matrix_q)
1484
1485 CALL timestop(handle)
1486
1487 END SUBROUTINE projector
1488
1489! **************************************************************************************************
1490!> \brief performs a tranformation of a matrix back to/into orthonormal basis
1491!> in case of P a scaling of 0.5 has to be applied for closed shell case
1492!> \param matrix matrix to be transformed
1493!> \param matrix_trafo transformation matrix
1494!> \param eps_filter filtering threshold for sparse matrices
1495!> \par History
1496!> 2012.05 created [Florian Schiffmann]
1497!> \author Florian Schiffmann
1498!>
1499! **************************************************************************************************
1500
1501 SUBROUTINE transform_m_orth(matrix, matrix_trafo, eps_filter)
1502 TYPE(dbcsr_type) :: matrix, matrix_trafo
1503 REAL(kind=dp) :: eps_filter
1504
1505 CHARACTER(LEN=*), PARAMETER :: routinen = 'transform_m_orth'
1506
1507 INTEGER :: handle
1508 TYPE(dbcsr_type) :: matrix_tmp, matrix_work
1509
1510 CALL timeset(routinen, handle)
1511
1512 CALL dbcsr_create(matrix_work, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1513 CALL dbcsr_create(matrix_tmp, template=matrix, matrix_type=dbcsr_type_no_symmetry)
1514
1515 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix, matrix_trafo, &
1516 0.0_dp, matrix_work, filter_eps=eps_filter)
1517 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_trafo, matrix_work, &
1518 0.0_dp, matrix_tmp, filter_eps=eps_filter)
1519 ! symmetrize results (this is again needed to make sure everything is stable)
1520 CALL dbcsr_transposed(matrix_work, matrix_tmp)
1521 CALL dbcsr_add(matrix_tmp, matrix_work, 0.5_dp, 0.5_dp)
1522 CALL dbcsr_copy(matrix, matrix_tmp)
1523
1524 ! Avoid the buildup of noisy blocks
1525 CALL dbcsr_filter(matrix, eps_filter)
1526
1527 CALL dbcsr_release(matrix_tmp)
1528 CALL dbcsr_release(matrix_work)
1529 CALL timestop(handle)
1530
1531 END SUBROUTINE transform_m_orth
1532
1533END MODULE ec_orth_solver
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:599
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_transposed
subroutine, public dbcsr_transposed(transposed, normal, shallow_data_copy, transpose_distribution, use_distribution)
...
Definition cp_dbcsr_api.F:1242
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_desymmetrize
subroutine, public dbcsr_desymmetrize(matrix_a, matrix_b)
...
Definition cp_dbcsr_api.F:525
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:658
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_checksum
real(kind=dp) function, public dbcsr_checksum(matrix, pos)
Calculates the checksum of a DBCSR matrix.
Definition cp_dbcsr_contrib.F:628
cp_dbcsr_contrib::dbcsr_add_on_diag
subroutine, public dbcsr_add_on_diag(matrix, alpha)
Adds the given scalar to the diagonal of the matrix. Reserves any missing diagonal blocks.
Definition cp_dbcsr_contrib.F:331
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_external_control
Routines to handle the external control of CP2K.
Definition cp_external_control.F:15
cp_external_control::external_control
subroutine, public external_control(should_stop, flag, globenv, target_time, start_time, force_check)
External manipulations during a run : when the <PROJECT_NAME>.EXIT_$runtype command is sent the progr...
Definition cp_external_control.F:91
ec_orth_solver
AO-based conjugate-gradient response solver routines.
Definition ec_orth_solver.F:15
ec_orth_solver::ec_response_ao
subroutine, public ec_response_ao(qs_env, p_env, matrix_hz, matrix_pz, matrix_wz, iounit, should_stop, silent)
AO-based conjugate gradient linear response solver. In goes the right hand side B of the equation AZ=...
Definition ec_orth_solver.F:369
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::precond_mlp
integer, parameter, public precond_mlp
Definition input_constants.F:533
input_constants::kg_tnadd_embed_ri
integer, parameter, public kg_tnadd_embed_ri
Definition input_constants.F:1204
input_constants::kg_tnadd_embed
integer, parameter, public kg_tnadd_embed
Definition input_constants.F:1204
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition input_constants.F:909
input_constants::ls_s_sqrt_proot
integer, parameter, public ls_s_sqrt_proot
Definition input_constants.F:1024
input_constants::ls_s_sqrt_ns
integer, parameter, public ls_s_sqrt_ns
Definition input_constants.F:1024
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
iterate_matrix
Routines useful for iterative matrix calculations.
Definition iterate_matrix.F:13
iterate_matrix::matrix_sqrt_newton_schulz
subroutine, public matrix_sqrt_newton_schulz(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged, iounit)
compute the sqrt of a matrix via the sign function and the corresponding Newton-Schulz iterations the...
Definition iterate_matrix.F:1627
iterate_matrix::matrix_sqrt_proot
subroutine, public matrix_sqrt_proot(matrix_sqrt, matrix_sqrt_inv, matrix, threshold, order, eps_lanczos, max_iter_lanczos, symmetrize, converged)
compute the sqrt of a matrix via the general algorithm for the p-th root of Richters et al....
Definition iterate_matrix.F:1880
kg_correction
Routines for a Kim-Gordon-like partitioning into molecular subunits.
Definition kg_correction.F:14
kg_correction::kg_ekin_subset
subroutine, public kg_ekin_subset(qs_env, ks_matrix, ekin_mol, calc_force, do_kernel, pmat_ext)
Calculates the subsystem Hohenberg-Kohn kinetic energy and the forces.
Definition kg_correction.F:91
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:124
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition machine.F:141
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::abnormal_value
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Definition mathlib.F:159
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_poisson_methods
Definition pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition pw_poisson_types.F:22
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:539
qs_fxc
Setup Routine for Fxc Potentials.
Definition qs_fxc.F:29
qs_fxc::qs_fxc_prep
subroutine, public qs_fxc_prep(qs_env, rho0_struct, xc_rho_set, xc_deriv_set, xc_section, pw_env_ext, is_triplet)
...
Definition qs_fxc.F:518
qs_fxc::qs_fxc_create
subroutine, public qs_fxc_create(qs_env, rho0_struct, rho1_struct, rho0_atom_set, xc_section, do_onecenter, fxc_rho, fxc_tau, rho1_atom_set, do_scale, is_triplet, spinflip, no_weights, uf_grid_results, pw_env_ext, kind_set_external, para_env_external, compute_virial, virial_xc)
...
Definition qs_fxc.F:118
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
qs_kpp1_env_types
basis types for the calculation of the perturbation of density theory.
Definition qs_kpp1_env_types.F:14
qs_linres_kernel
linres kernel functions
Definition qs_linres_kernel.F:14
qs_linres_kernel::apply_xc_admm
subroutine, public apply_xc_admm(qs_env, p_env)
...
Definition qs_linres_kernel.F:959
qs_linres_kernel::apply_hfx
subroutine, public apply_hfx(qs_env, p_env)
Update action of TDDFPT operator on trial vectors by adding exact-exchange term.
Definition qs_linres_kernel.F:795
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_p_env_methods
Utility functions for the perturbation calculations.
Definition qs_p_env_methods.F:16
qs_p_env_methods::p_env_finish_kpp1
subroutine, public p_env_finish_kpp1(qs_env, p_env)
...
Definition qs_p_env_methods.F:853
qs_p_env_methods::p_env_update_rho
subroutine, public p_env_update_rho(p_env, qs_env)
...
Definition qs_p_env_methods.F:395
qs_p_env_methods::p_env_check_i_alloc
subroutine, public p_env_check_i_alloc(p_env, qs_env)
checks that the intenal storage is allocated, and allocs it if needed
Definition qs_p_env_methods.F:314
qs_p_env_types
basis types for the calculation of the perturbation of density theory.
Definition qs_p_env_types.F:14
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
cp_control_types::dft_control_type
Definition cp_control_types.F:758
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_p_type
to create arrays of pools
Definition pw_pool_types.F:135
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kpp1_env_types::qs_kpp1_env_type
environment that keeps the informations and temporary val to build the kpp1 kernel matrix
Definition qs_kpp1_env_types.F:40
qs_linres_types::linres_control_type
General settings for linear response calculations.
Definition qs_linres_types.F:53
qs_p_env_types::qs_p_env_type
Represent a qs system that is perturbed. Can calculate the linear operator and the rhs of the system ...
Definition qs_p_env_types.F:58
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62