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