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