(git:b279b6b)
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 ! **************************************************************************************************
15 MODULE 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 
85 CONTAINS
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 
1517 END 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
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