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