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