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