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