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 Routines to calculate CPHF like update and solve Z-vector equation
10 : !> for MP2 gradients (only GPW)
11 : !> \par History
12 : !> 11.2013 created [Mauro Del Ben]
13 : ! **************************************************************************************************
14 : MODULE mp2_cphf
15 : USE admm_methods, ONLY: admm_projection_derivative
16 : USE admm_types, ONLY: admm_type,&
17 : get_admm_env
18 : USE cell_types, ONLY: cell_type
19 : USE cp_blacs_env, ONLY: cp_blacs_env_type
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
22 : dbcsr_copy,&
23 : dbcsr_p_type,&
24 : dbcsr_release,&
25 : dbcsr_scale,&
26 : dbcsr_set
27 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
28 : cp_dbcsr_plus_fm_fm_t,&
29 : dbcsr_allocate_matrix_set,&
30 : dbcsr_deallocate_matrix_set
31 : USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_to_full
32 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
33 : cp_fm_struct_p_type,&
34 : cp_fm_struct_release,&
35 : cp_fm_struct_type
36 : USE cp_fm_types, ONLY: cp_fm_create,&
37 : cp_fm_get_info,&
38 : cp_fm_release,&
39 : cp_fm_set_all,&
40 : cp_fm_to_fm_submat,&
41 : cp_fm_type
42 : USE cp_log_handling, ONLY: cp_get_default_logger,&
43 : cp_logger_type
44 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
45 : cp_print_key_unit_nr
46 : USE hfx_admm_utils, ONLY: tddft_hfx_matrix
47 : USE hfx_derivatives, ONLY: derivatives_four_center
48 : USE hfx_exx, ONLY: add_exx_to_rhs
49 : USE hfx_ri, ONLY: hfx_ri_update_forces
50 : USE hfx_types, ONLY: alloc_containers,&
51 : hfx_container_type,&
52 : hfx_init_container,&
53 : hfx_type
54 : USE input_constants, ONLY: do_admm_aux_exch_func_none,&
55 : ot_precond_full_all,&
56 : z_solver_cg,&
57 : z_solver_pople,&
58 : z_solver_richardson,&
59 : z_solver_sd
60 : USE input_section_types, ONLY: section_vals_get,&
61 : section_vals_get_subs_vals,&
62 : section_vals_type
63 : USE kahan_sum, ONLY: accurate_dot_product
64 : USE kinds, ONLY: dp
65 : USE linear_systems, ONLY: solve_system
66 : USE machine, ONLY: m_flush,&
67 : m_walltime
68 : USE mathconstants, ONLY: fourpi
69 : USE message_passing, ONLY: mp_para_env_type
70 : USE mp2_types, ONLY: mp2_type,&
71 : ri_rpa_method_gpw
72 : USE parallel_gemm_api, ONLY: parallel_gemm
73 : USE pw_env_types, ONLY: pw_env_get,&
74 : pw_env_type
75 : USE pw_methods, ONLY: pw_axpy,&
76 : pw_copy,&
77 : pw_derive,&
78 : pw_integral_ab,&
79 : pw_scale,&
80 : pw_transfer,&
81 : pw_zero
82 : USE pw_poisson_methods, ONLY: pw_poisson_solve
83 : USE pw_poisson_types, ONLY: pw_poisson_type
84 : USE pw_pool_types, ONLY: pw_pool_type
85 : USE pw_types, ONLY: pw_c1d_gs_type,&
86 : pw_r3d_rs_type
87 : USE qs_2nd_kernel_ao, ONLY: apply_2nd_order_kernel
88 : USE qs_core_matrices, ONLY: core_matrices,&
89 : kinetic_energy_matrix
90 : USE qs_density_matrices, ONLY: calculate_whz_matrix
91 : USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
92 : USE qs_dispersion_types, ONLY: qs_dispersion_type
93 : USE qs_energy_types, ONLY: qs_energy_type
94 : USE qs_environment_types, ONLY: get_qs_env,&
95 : qs_environment_type,&
96 : set_qs_env
97 : USE qs_force_types, ONLY: deallocate_qs_force,&
98 : qs_force_type,&
99 : sum_qs_force,&
100 : zero_qs_force
101 : USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
102 : integrate_v_rspace
103 : USE qs_ks_reference, ONLY: ks_ref_potential
104 : USE qs_ks_types, ONLY: qs_ks_env_type,&
105 : set_ks_env
106 : USE qs_linres_types, ONLY: linres_control_type
107 : USE qs_mo_types, ONLY: get_mo_set,&
108 : mo_set_type
109 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
110 : USE qs_overlap, ONLY: build_overlap_matrix
111 : USE qs_p_env_methods, ONLY: p_env_check_i_alloc,&
112 : p_env_create,&
113 : p_env_psi0_changed,&
114 : p_env_update_rho
115 : USE qs_p_env_types, ONLY: p_env_release,&
116 : qs_p_env_type
117 : USE qs_rho_types, ONLY: qs_rho_get,&
118 : qs_rho_type
119 : USE task_list_types, ONLY: task_list_type
120 : USE virial_types, ONLY: virial_type,&
121 : zero_virial
122 :
123 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
124 :
125 : #include "./base/base_uses.f90"
126 :
127 : IMPLICIT NONE
128 :
129 : PRIVATE
130 :
131 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_cphf'
132 : LOGICAL, PARAMETER, PRIVATE :: debug_forces = .TRUE.
133 :
134 : PUBLIC :: solve_z_vector_eq, update_mp2_forces
135 :
136 : CONTAINS
137 :
138 : ! **************************************************************************************************
139 : !> \brief Solve Z-vector equations necessary for the calculation of the MP2
140 : !> gradients, in order to be consistent here the parameters for the
141 : !> calculation of the CPHF like updats have to be exactly equal to the
142 : !> SCF case
143 : !> \param qs_env ...
144 : !> \param mp2_env ...
145 : !> \param para_env ...
146 : !> \param dft_control ...
147 : !> \param mo_coeff ...
148 : !> \param homo ...
149 : !> \param Eigenval ...
150 : !> \param unit_nr ...
151 : !> \author Mauro Del Ben, Vladimir Rybkin
152 : ! **************************************************************************************************
153 272 : SUBROUTINE solve_z_vector_eq(qs_env, mp2_env, para_env, dft_control, &
154 272 : mo_coeff, homo, Eigenval, unit_nr)
155 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
156 : TYPE(mp2_type), INTENT(INOUT) :: mp2_env
157 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
158 : TYPE(dft_control_type), INTENT(IN), POINTER :: dft_control
159 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff
160 : INTEGER, DIMENSION(:), INTENT(IN) :: homo
161 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
162 : INTEGER, INTENT(IN) :: unit_nr
163 :
164 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_z_vector_eq'
165 :
166 : INTEGER :: bin, handle, handle2, i, i_global, i_thread, iiB, irep, ispin, j_global, jjB, &
167 : my_bin_size, n_rep_hf, n_threads, nao, ncol_local, nmo, nrow_local, nspins, transf_type_in
168 272 : INTEGER, ALLOCATABLE, DIMENSION(:) :: virtual
169 272 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
170 : LOGICAL :: alpha_beta, do_dynamic_load_balancing, &
171 : do_exx, do_hfx, restore_p_screen
172 : REAL(KIND=dp) :: focc
173 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
174 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
175 : TYPE(cp_fm_type) :: fm_back, fm_G_mu_nu, fm_mo_mo
176 272 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: L_jb, mo_coeff_o, mo_coeff_v, P_ia, &
177 272 : P_mo, W_Mo
178 272 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p_mp2, &
179 272 : matrix_p_mp2_admm, matrix_s, P_mu_nu, &
180 272 : rho1_ao, rho_ao, rho_ao_aux_fit
181 272 : TYPE(hfx_container_type), DIMENSION(:), POINTER :: integral_containers
182 : TYPE(hfx_container_type), POINTER :: maxval_container
183 : TYPE(hfx_type), POINTER :: actual_x_data
184 : TYPE(linres_control_type), POINTER :: linres_control
185 : TYPE(qs_ks_env_type), POINTER :: ks_env
186 : TYPE(qs_p_env_type), POINTER :: p_env
187 : TYPE(qs_rho_type), POINTER :: rho, rho_aux_fit
188 : TYPE(section_vals_type), POINTER :: hfx_section, hfx_sections, input
189 :
190 272 : CALL timeset(routineN, handle)
191 :
192 : ! start collecting stuff
193 272 : CALL cp_fm_get_info(mo_coeff(1), nrow_global=nao, ncol_global=nmo)
194 272 : CPASSERT(SIZE(eigenval, 1) == nmo)
195 272 : CPASSERT(nao >= nmo)
196 272 : NULLIFY (input, matrix_s, blacs_env, rho, &
197 272 : matrix_p_mp2, matrix_p_mp2_admm, matrix_ks)
198 : CALL get_qs_env(qs_env, &
199 : ks_env=ks_env, &
200 : input=input, &
201 : matrix_s=matrix_s, &
202 : matrix_ks=matrix_ks, &
203 : matrix_p_mp2=matrix_p_mp2, &
204 : matrix_p_mp2_admm=matrix_p_mp2_admm, &
205 : blacs_env=blacs_env, &
206 272 : rho=rho)
207 :
208 272 : CALL qs_rho_get(rho, rho_ao=rho_ao)
209 :
210 : ! Get number of relevant spin states
211 272 : nspins = dft_control%nspins
212 272 : alpha_beta = (nspins == 2)
213 :
214 272 : CALL MOVE_ALLOC(mp2_env%ri_grad%P_mo, P_mo)
215 272 : CALL MOVE_ALLOC(mp2_env%ri_grad%W_mo, W_mo)
216 272 : CALL MOVE_ALLOC(mp2_env%ri_grad%L_jb, L_jb)
217 :
218 816 : ALLOCATE (virtual(nspins))
219 630 : virtual(:) = nmo - homo(:)
220 :
221 272 : NULLIFY (P_mu_nu)
222 272 : CALL dbcsr_allocate_matrix_set(P_mu_nu, nspins)
223 630 : DO ispin = 1, nspins
224 358 : ALLOCATE (P_mu_nu(ispin)%matrix)
225 358 : CALL dbcsr_copy(P_mu_nu(ispin)%matrix, rho_ao(1)%matrix, name="P_mu_nu")
226 630 : CALL dbcsr_set(P_mu_nu(ispin)%matrix, 0.0_dp)
227 : END DO
228 :
229 272 : NULLIFY (fm_struct_tmp)
230 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
231 272 : nrow_global=nao, ncol_global=nao)
232 272 : CALL cp_fm_create(fm_G_mu_nu, fm_struct_tmp, name="G_mu_nu")
233 272 : CALL cp_fm_create(fm_back, fm_struct_tmp, name="fm_back")
234 272 : CALL cp_fm_struct_release(fm_struct_tmp)
235 272 : CALL cp_fm_set_all(fm_G_mu_nu, 0.0_dp)
236 272 : CALL cp_fm_set_all(fm_back, 0.0_dp)
237 :
238 1804 : ALLOCATE (mo_coeff_o(nspins), mo_coeff_v(nspins))
239 630 : DO ispin = 1, nspins
240 358 : NULLIFY (fm_struct_tmp)
241 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
242 358 : nrow_global=nao, ncol_global=homo(ispin))
243 358 : CALL cp_fm_create(mo_coeff_o(ispin), fm_struct_tmp, name="mo_coeff_o")
244 358 : CALL cp_fm_struct_release(fm_struct_tmp)
245 358 : CALL cp_fm_set_all(mo_coeff_o(ispin), 0.0_dp)
246 : CALL cp_fm_to_fm_submat(msource=mo_coeff(ispin), mtarget=mo_coeff_o(ispin), &
247 : nrow=nao, ncol=homo(ispin), &
248 : s_firstrow=1, s_firstcol=1, &
249 358 : t_firstrow=1, t_firstcol=1)
250 :
251 358 : NULLIFY (fm_struct_tmp)
252 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
253 358 : nrow_global=nao, ncol_global=virtual(ispin))
254 358 : CALL cp_fm_create(mo_coeff_v(ispin), fm_struct_tmp, name="mo_coeff_v")
255 358 : CALL cp_fm_struct_release(fm_struct_tmp)
256 358 : CALL cp_fm_set_all(mo_coeff_v(ispin), 0.0_dp)
257 : CALL cp_fm_to_fm_submat(msource=mo_coeff(ispin), mtarget=mo_coeff_v(ispin), &
258 : nrow=nao, ncol=virtual(ispin), &
259 : s_firstrow=1, s_firstcol=homo(ispin) + 1, &
260 630 : t_firstrow=1, t_firstcol=1)
261 : END DO
262 :
263 : ! hfx section
264 272 : NULLIFY (hfx_sections)
265 272 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
266 272 : CALL section_vals_get(hfx_sections, explicit=do_hfx, n_repetition=n_rep_hf)
267 272 : IF (do_hfx) THEN
268 : ! here we check if we have to reallocate the HFX container
269 192 : IF (mp2_env%ri_grad%free_hfx_buffer .AND. (.NOT. qs_env%x_data(1, 1)%do_hfx_ri)) THEN
270 102 : CALL timeset(routineN//"_alloc_hfx", handle2)
271 102 : n_threads = 1
272 102 : !$ n_threads = omp_get_max_threads()
273 :
274 204 : DO irep = 1, n_rep_hf
275 306 : DO i_thread = 0, n_threads - 1
276 102 : actual_x_data => qs_env%x_data(irep, i_thread + 1)
277 :
278 102 : do_dynamic_load_balancing = .TRUE.
279 102 : IF (n_threads == 1 .OR. actual_x_data%memory_parameter%do_disk_storage) do_dynamic_load_balancing = .FALSE.
280 :
281 : IF (do_dynamic_load_balancing) THEN
282 0 : my_bin_size = SIZE(actual_x_data%distribution_energy)
283 : ELSE
284 102 : my_bin_size = 1
285 : END IF
286 :
287 204 : IF (.NOT. actual_x_data%memory_parameter%do_all_on_the_fly) THEN
288 102 : CALL alloc_containers(actual_x_data%store_ints, my_bin_size)
289 :
290 204 : DO bin = 1, my_bin_size
291 102 : maxval_container => actual_x_data%store_ints%maxval_container(bin)
292 102 : integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
293 102 : CALL hfx_init_container(maxval_container, actual_x_data%memory_parameter%actual_memory_usage, .FALSE.)
294 6732 : DO i = 1, 64
295 : CALL hfx_init_container(integral_containers(i), &
296 6630 : actual_x_data%memory_parameter%actual_memory_usage, .FALSE.)
297 : END DO
298 : END DO
299 : END IF
300 : END DO
301 : END DO
302 102 : CALL timestop(handle2)
303 : END IF
304 :
305 : ! set up parameters for P_screening
306 192 : restore_p_screen = qs_env%x_data(1, 1)%screening_parameter%do_initial_p_screening
307 192 : IF (qs_env%x_data(1, 1)%screening_parameter%do_initial_p_screening) THEN
308 8 : IF (mp2_env%ri_grad%free_hfx_buffer) THEN
309 8 : mp2_env%p_screen = .FALSE.
310 : ELSE
311 0 : mp2_env%p_screen = .TRUE.
312 : END IF
313 : END IF
314 : END IF
315 :
316 : ! Add exx part for RPA
317 272 : do_exx = .FALSE.
318 272 : IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
319 24 : hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
320 24 : CALL section_vals_get(hfx_section, explicit=do_exx)
321 : END IF
322 272 : IF (do_exx) THEN
323 : CALL add_exx_to_rhs(rhs=P_mu_nu, &
324 : qs_env=qs_env, &
325 : ext_hfx_section=hfx_section, &
326 : x_data=mp2_env%ri_rpa%x_data, &
327 : recalc_integrals=.FALSE., &
328 : do_admm=mp2_env%ri_rpa%do_admm, &
329 : do_exx=do_exx, &
330 8 : reuse_hfx=mp2_env%ri_rpa%reuse_hfx)
331 :
332 8 : focc = 1.0_dp
333 8 : IF (nspins == 1) focc = 2.0_dp
334 : !focc = 0.0_dp
335 16 : DO ispin = 1, nspins
336 8 : CALL dbcsr_add(P_mu_nu(ispin)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, -1.0_dp)
337 8 : CALL copy_dbcsr_to_fm(matrix=P_mu_nu(ispin)%matrix, fm=fm_G_mu_nu)
338 : CALL parallel_gemm("N", "N", nao, homo(ispin), nao, 1.0_dp, &
339 8 : fm_G_mu_nu, mo_coeff_o(ispin), 0.0_dp, fm_back)
340 : CALL parallel_gemm("T", "N", homo(ispin), virtual(ispin), nao, focc, &
341 8 : fm_back, mo_coeff_v(ispin), 1.0_dp, L_jb(ispin))
342 : CALL parallel_gemm("T", "N", homo(ispin), homo(ispin), nao, -focc, &
343 16 : fm_back, mo_coeff_o(ispin), 1.0_dp, W_mo(ispin))
344 : END DO
345 : END IF
346 :
347 : ! Prepare arrays for linres code
348 : NULLIFY (linres_control)
349 272 : ALLOCATE (linres_control)
350 272 : linres_control%do_kernel = .TRUE.
351 : linres_control%lr_triplet = .FALSE.
352 : linres_control%linres_restart = .FALSE.
353 272 : linres_control%max_iter = mp2_env%ri_grad%cphf_max_num_iter
354 272 : linres_control%eps = mp2_env%ri_grad%cphf_eps_conv
355 272 : linres_control%eps_filter = mp2_env%mp2_gpw%eps_filter
356 272 : linres_control%restart_every = 50
357 272 : linres_control%preconditioner_type = ot_precond_full_all
358 272 : linres_control%energy_gap = 0.02_dp
359 :
360 : NULLIFY (p_env)
361 1904 : ALLOCATE (p_env)
362 272 : CALL p_env_create(p_env, qs_env, p1_option=P_mu_nu, orthogonal_orbitals=.TRUE., linres_control=linres_control)
363 272 : CALL set_qs_env(qs_env, linres_control=linres_control)
364 272 : CALL p_env_psi0_changed(p_env, qs_env)
365 272 : p_env%new_preconditioner = .TRUE.
366 272 : CALL p_env_check_i_alloc(p_env, qs_env)
367 272 : mp2_env%ri_grad%p_env => p_env
368 :
369 : ! update Lagrangian with the CPHF like update, occ-occ block, first call (recompute hfx integrals if needed)
370 272 : transf_type_in = 1
371 : ! In alpha-beta case, L_bj_alpha has Coulomb and XC alpha-alpha part
372 : ! and (only) Coulomb alpha-beta part and vice versa.
373 :
374 : ! Complete in closed shell case, alpha-alpha (Coulomb and XC)
375 : ! part of L_bj(alpha) for open shell
376 :
377 : CALL cphf_like_update(qs_env, mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
378 : P_mo, fm_G_mu_nu, fm_back, transf_type_in, L_jb, &
379 : recalc_hfx_integrals=(.NOT. do_exx .AND. mp2_env%ri_grad%free_hfx_buffer) &
380 364 : .OR. (do_exx .AND. .NOT. mp2_env%ri_rpa%reuse_hfx))
381 :
382 : ! at this point Lagrangian is completed ready to solve the Z-vector equations
383 : ! P_ia will contain the solution of these equations
384 902 : ALLOCATE (P_ia(nspins))
385 630 : DO ispin = 1, nspins
386 358 : NULLIFY (fm_struct_tmp)
387 : CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
388 358 : nrow_global=homo(ispin), ncol_global=virtual(ispin))
389 358 : CALL cp_fm_create(P_ia(ispin), fm_struct_tmp, name="P_ia")
390 358 : CALL cp_fm_struct_release(fm_struct_tmp)
391 630 : CALL cp_fm_set_all(P_ia(ispin), 0.0_dp)
392 : END DO
393 :
394 : CALL solve_z_vector_eq_low(qs_env, mp2_env, unit_nr, &
395 : mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
396 272 : L_jb, fm_G_mu_nu, fm_back, P_ia)
397 :
398 : ! release fm stuff
399 272 : CALL cp_fm_release(fm_G_mu_nu)
400 272 : CALL cp_fm_release(L_jb)
401 272 : CALL cp_fm_release(mo_coeff_o)
402 272 : CALL cp_fm_release(mo_coeff_v)
403 :
404 272 : CALL cp_fm_create(fm_mo_mo, P_mo(1)%matrix_struct)
405 :
406 630 : DO ispin = 1, nspins
407 : ! update the MP2-MO density matrix with the occ-virt block
408 : CALL cp_fm_to_fm_submat(msource=P_ia(ispin), mtarget=P_mo(ispin), &
409 : nrow=homo(ispin), ncol=virtual(ispin), &
410 : s_firstrow=1, s_firstcol=1, &
411 358 : t_firstrow=1, t_firstcol=homo(ispin) + 1)
412 : ! transpose P_MO matrix (easy way to symmetrize)
413 358 : CALL cp_fm_set_all(fm_mo_mo, 0.0_dp)
414 : ! P_mo now is ready
415 630 : CALL cp_fm_uplo_to_full(matrix=P_mo(ispin), work=fm_mo_mo)
416 : END DO
417 272 : CALL cp_fm_release(P_ia)
418 272 : CALL cp_fm_release(fm_mo_mo)
419 :
420 : ! do the final update to the energy weighted matrix W_MO
421 630 : DO ispin = 1, nspins
422 : CALL cp_fm_get_info(matrix=W_mo(ispin), &
423 : nrow_local=nrow_local, &
424 : ncol_local=ncol_local, &
425 : row_indices=row_indices, &
426 358 : col_indices=col_indices)
427 7402 : DO jjB = 1, ncol_local
428 6772 : j_global = col_indices(jjB)
429 7130 : IF (j_global <= homo(ispin)) THEN
430 13994 : DO iiB = 1, nrow_local
431 12692 : i_global = row_indices(iiB)
432 : W_mo(ispin)%local_data(iiB, jjB) = W_mo(ispin)%local_data(iiB, jjB) &
433 12692 : - P_mo(ispin)%local_data(iiB, jjB)*Eigenval(j_global, ispin)
434 12692 : IF (i_global == j_global .AND. nspins == 1) W_mo(ispin)%local_data(iiB, jjB) = &
435 359 : W_mo(ispin)%local_data(iiB, jjB) - 2.0_dp*Eigenval(j_global, ispin)
436 12692 : IF (i_global == j_global .AND. nspins == 2) W_mo(ispin)%local_data(iiB, jjB) = &
437 1594 : W_mo(ispin)%local_data(iiB, jjB) - Eigenval(j_global, ispin)
438 : END DO
439 : ELSE
440 67764 : DO iiB = 1, nrow_local
441 62294 : i_global = row_indices(iiB)
442 67764 : IF (i_global <= homo(ispin)) THEN
443 : ! virt-occ
444 : W_mo(ispin)%local_data(iiB, jjB) = W_mo(ispin)%local_data(iiB, jjB) &
445 10193 : - P_mo(ispin)%local_data(iiB, jjB)*Eigenval(i_global, ispin)
446 : ELSE
447 : ! virt-virt
448 : W_mo(ispin)%local_data(iiB, jjB) = W_mo(ispin)%local_data(iiB, jjB) &
449 52101 : - P_mo(ispin)%local_data(iiB, jjB)*Eigenval(j_global, ispin)
450 : END IF
451 : END DO
452 : END IF
453 : END DO
454 : END DO
455 :
456 : ! create the MP2 energy weighted density matrix
457 272 : NULLIFY (p_env%w1)
458 272 : CALL dbcsr_allocate_matrix_set(p_env%w1, 1)
459 272 : ALLOCATE (p_env%w1(1)%matrix)
460 : CALL dbcsr_copy(p_env%w1(1)%matrix, matrix_s(1)%matrix, &
461 272 : name="W MATRIX MP2")
462 272 : CALL dbcsr_set(p_env%w1(1)%matrix, 0.0_dp)
463 :
464 : ! backtnsform the collected parts of the energy-weighted density matrix into AO basis
465 630 : DO ispin = 1, nspins
466 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, 1.0_dp, &
467 358 : mo_coeff(ispin), W_mo(ispin), 0.0_dp, fm_back)
468 630 : CALL cp_dbcsr_plus_fm_fm_t(p_env%w1(1)%matrix, fm_back, mo_coeff(ispin), nmo, 1.0_dp, .TRUE., 1)
469 : END DO
470 272 : CALL cp_fm_release(W_mo)
471 :
472 272 : CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
473 :
474 630 : DO ispin = 1, nspins
475 358 : CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
476 :
477 : CALL parallel_gemm('N', 'N', nao, nmo, nmo, 1.0_dp, &
478 358 : mo_coeff(ispin), P_mo(ispin), 0.0_dp, fm_back)
479 358 : CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff(ispin), nmo, 1.0_dp, .TRUE.)
480 :
481 630 : CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix)
482 : END DO
483 272 : CALL cp_fm_release(P_mo)
484 272 : CALL cp_fm_release(fm_back)
485 :
486 272 : CALL p_env_update_rho(p_env, qs_env)
487 :
488 : ! create mp2 DBCSR density
489 272 : CALL dbcsr_allocate_matrix_set(matrix_p_mp2, nspins)
490 630 : DO ispin = 1, nspins
491 358 : ALLOCATE (matrix_p_mp2(ispin)%matrix)
492 : CALL dbcsr_copy(matrix_p_mp2(ispin)%matrix, p_env%p1(ispin)%matrix, &
493 630 : name="P MATRIX MP2")
494 : END DO
495 :
496 272 : IF (dft_control%do_admm) THEN
497 40 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux_fit)
498 40 : CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
499 :
500 : ! create mp2 DBCSR density in auxiliary basis
501 40 : CALL dbcsr_allocate_matrix_set(matrix_p_mp2_admm, nspins)
502 98 : DO ispin = 1, nspins
503 58 : ALLOCATE (matrix_p_mp2_admm(ispin)%matrix)
504 : CALL dbcsr_copy(matrix_p_mp2_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix, &
505 98 : name="P MATRIX MP2 ADMM")
506 : END DO
507 : END IF
508 :
509 272 : CALL set_ks_env(ks_env, matrix_p_mp2=matrix_p_mp2, matrix_p_mp2_admm=matrix_p_mp2_admm)
510 :
511 : ! We will need one more hfx calculation for HF gradient part
512 272 : mp2_env%not_last_hfx = .FALSE.
513 272 : mp2_env%p_screen = restore_p_screen
514 :
515 272 : CALL timestop(handle)
516 :
517 1360 : END SUBROUTINE solve_z_vector_eq
518 :
519 : ! **************************************************************************************************
520 : !> \brief Here we performe the CPHF like update using GPW,
521 : !> transf_type_in defines the type of transformation for the matrix in input
522 : !> transf_type_in = 1 -> occ-occ back transformation
523 : !> transf_type_in = 2 -> virt-virt back transformation
524 : !> transf_type_in = 3 -> occ-virt back transformation including the
525 : !> eigenvalues energy differences for the diagonal elements
526 : !> \param qs_env ...
527 : !> \param mo_coeff_o ...
528 : !> \param mo_coeff_v ...
529 : !> \param Eigenval ...
530 : !> \param p_env ...
531 : !> \param fm_mo ...
532 : !> \param fm_ao ...
533 : !> \param fm_back ...
534 : !> \param transf_type_in ...
535 : !> \param fm_mo_out ...
536 : !> \param recalc_hfx_integrals ...
537 : !> \author Mauro Del Ben, Vladimir Rybkin
538 : ! **************************************************************************************************
539 1512 : SUBROUTINE cphf_like_update(qs_env, mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
540 1512 : fm_mo, fm_ao, fm_back, transf_type_in, &
541 1512 : fm_mo_out, recalc_hfx_integrals)
542 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
543 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff_o, mo_coeff_v
544 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
545 : TYPE(qs_p_env_type) :: p_env
546 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mo
547 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_ao
548 : TYPE(cp_fm_type), INTENT(IN) :: fm_back
549 : INTEGER, INTENT(IN) :: transf_type_in
550 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: fm_mo_out
551 : LOGICAL, INTENT(IN), OPTIONAL :: recalc_hfx_integrals
552 :
553 : CHARACTER(LEN=*), PARAMETER :: routineN = 'cphf_like_update'
554 :
555 : INTEGER :: handle, homo, i_global, iiB, ispin, &
556 : j_global, jjB, nao, ncol_local, &
557 : nrow_local, nspins, virtual
558 1512 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
559 1512 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho1_ao
560 :
561 1512 : CALL timeset(routineN, handle)
562 :
563 1512 : nspins = SIZE(Eigenval, 2)
564 :
565 1512 : CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
566 :
567 : ! Determine the first-order density matrices in AO basis
568 3488 : DO ispin = 1, nspins
569 1976 : CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
570 :
571 1976 : CALL cp_fm_get_info(mo_coeff_o(ispin), nrow_global=nao, ncol_global=homo)
572 1976 : CALL cp_fm_get_info(mo_coeff_v(ispin), nrow_global=nao, ncol_global=virtual)
573 :
574 : ASSOCIATE (mat_in => fm_mo(ispin))
575 :
576 : ! perform back transformation
577 2334 : SELECT CASE (transf_type_in)
578 : CASE (1)
579 : ! occ-occ block
580 : CALL parallel_gemm('N', 'N', nao, homo, homo, 1.0_dp, &
581 358 : mo_coeff_o(ispin), mat_in, 0.0_dp, fm_back)
582 358 : CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_o(ispin), homo, 1.0_dp, .TRUE.)
583 : ! virt-virt block
584 : CALL parallel_gemm('N', 'N', nao, virtual, virtual, 1.0_dp, &
585 : mo_coeff_v(ispin), mat_in, 0.0_dp, fm_back, &
586 : b_first_col=homo + 1, &
587 358 : b_first_row=homo + 1)
588 358 : CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_v(ispin), virtual, 1.0_dp, .TRUE.)
589 :
590 : CASE (3)
591 : ! virt-occ blocks
592 : CALL parallel_gemm('N', 'N', nao, virtual, homo, 1.0_dp, &
593 1618 : mo_coeff_o(ispin), mat_in, 0.0_dp, fm_back)
594 1618 : CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_v(ispin), virtual, 1.0_dp, .TRUE.)
595 : ! and symmetrize (here again multiply instead of transposing)
596 : CALL parallel_gemm('N', 'T', nao, homo, virtual, 1.0_dp, &
597 1618 : mo_coeff_v(ispin), mat_in, 0.0_dp, fm_back)
598 3594 : CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_o(ispin), homo, 1.0_dp, .TRUE.)
599 :
600 : CASE DEFAULT
601 : ! nothing
602 : END SELECT
603 : END ASSOCIATE
604 :
605 5464 : CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix)
606 : END DO
607 :
608 1512 : CALL p_env_update_rho(p_env, qs_env)
609 :
610 1512 : CALL apply_2nd_order_kernel(qs_env, p_env, recalc_hfx_integrals)
611 :
612 3488 : DO ispin = 1, nspins
613 1976 : CALL cp_fm_get_info(mo_coeff_o(ispin), nrow_global=nao, ncol_global=homo)
614 1976 : CALL cp_fm_get_info(mo_coeff_v(ispin), nrow_global=nao, ncol_global=virtual)
615 :
616 1976 : IF (transf_type_in == 3) THEN
617 :
618 : ! scale for the orbital energy differences for the diagonal elements
619 : CALL cp_fm_get_info(matrix=fm_mo_out(ispin), &
620 : nrow_local=nrow_local, &
621 : ncol_local=ncol_local, &
622 : row_indices=row_indices, &
623 1618 : col_indices=col_indices)
624 28110 : DO jjB = 1, ncol_local
625 26492 : j_global = col_indices(jjB)
626 79179 : DO iiB = 1, nrow_local
627 51069 : i_global = row_indices(iiB)
628 : fm_mo_out(ispin)%local_data(iiB, jjB) = fm_mo(ispin)%local_data(iiB, jjB)* &
629 77561 : (Eigenval(j_global + homo, ispin) - Eigenval(i_global, ispin))
630 : END DO
631 : END DO
632 : END IF
633 :
634 : ! copy back to fm
635 1976 : CALL cp_fm_set_all(fm_ao, 0.0_dp)
636 1976 : CALL copy_dbcsr_to_fm(matrix=p_env%kpp1(ispin)%matrix, fm=fm_ao)
637 1976 : CALL cp_fm_set_all(fm_back, 0.0_dp)
638 1976 : CALL cp_fm_uplo_to_full(fm_ao, fm_back)
639 :
640 3488 : ASSOCIATE (mat_out => fm_mo_out(ispin))
641 :
642 : ! transform to MO basis, here we always sum the result into the input matrix
643 :
644 : ! occ-virt block
645 : CALL parallel_gemm('T', 'N', homo, nao, nao, 1.0_dp, &
646 1976 : mo_coeff_o(ispin), fm_ao, 0.0_dp, fm_back)
647 : CALL parallel_gemm('N', 'N', homo, virtual, nao, 1.0_dp, &
648 1976 : fm_back, mo_coeff_v(ispin), 1.0_dp, mat_out)
649 : END ASSOCIATE
650 : END DO
651 :
652 1512 : CALL timestop(handle)
653 :
654 1512 : END SUBROUTINE cphf_like_update
655 :
656 : ! **************************************************************************************************
657 : !> \brief Low level subroutine for the iterative solution of a large
658 : !> system of linear equation
659 : !> \param qs_env ...
660 : !> \param mp2_env ...
661 : !> \param unit_nr ...
662 : !> \param mo_coeff_o ...
663 : !> \param mo_coeff_v ...
664 : !> \param Eigenval ...
665 : !> \param p_env ...
666 : !> \param L_jb ...
667 : !> \param fm_G_mu_nu ...
668 : !> \param fm_back ...
669 : !> \param P_ia ...
670 : !> \author Mauro Del Ben, Vladimir Rybkin
671 : ! **************************************************************************************************
672 272 : SUBROUTINE solve_z_vector_eq_low(qs_env, mp2_env, unit_nr, &
673 544 : mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
674 272 : L_jb, fm_G_mu_nu, fm_back, P_ia)
675 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
676 : TYPE(mp2_type), INTENT(IN) :: mp2_env
677 : INTEGER, INTENT(IN) :: unit_nr
678 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff_o, mo_coeff_v
679 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
680 : TYPE(qs_p_env_type), INTENT(IN), POINTER :: p_env
681 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: L_jb
682 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_G_mu_nu
683 : TYPE(cp_fm_type), INTENT(IN) :: fm_back
684 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: P_ia
685 :
686 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_z_vector_eq_low'
687 :
688 : INTEGER :: handle, i_global, iiB, ispin, j_global, &
689 : jjB, ncol_local, nmo, nrow_local, &
690 : nspins, virtual
691 272 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
692 272 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: precond
693 :
694 272 : CALL timeset(routineN, handle)
695 :
696 272 : nmo = SIZE(eigenval, 1)
697 272 : nspins = SIZE(eigenval, 2)
698 :
699 : ! Pople method
700 : ! change sign to L_jb
701 630 : DO ispin = 1, nspins
702 16320 : L_jb(ispin)%local_data(:, :) = -L_jb(ispin)%local_data(:, :)
703 : END DO
704 :
705 : ! create fm structure
706 1174 : ALLOCATE (precond(nspins))
707 630 : DO ispin = 1, nspins
708 : ! create preconditioner (for now only orbital energy differences)
709 358 : CALL cp_fm_create(precond(ispin), P_ia(ispin)%matrix_struct, name="precond")
710 358 : CALL cp_fm_set_all(precond(ispin), 1.0_dp)
711 : CALL cp_fm_get_info(matrix=precond(ispin), &
712 : nrow_local=nrow_local, &
713 : ncol_local=ncol_local, &
714 : row_indices=row_indices, &
715 : col_indices=col_indices, &
716 358 : ncol_global=virtual)
717 6100 : DO jjB = 1, ncol_local
718 5470 : j_global = col_indices(jjB)
719 16021 : DO iiB = 1, nrow_local
720 10193 : i_global = row_indices(iiB)
721 15663 : precond(ispin)%local_data(iiB, jjB) = 1.0_dp/(Eigenval(j_global + nmo - virtual, ispin) - Eigenval(i_global, ispin))
722 : END DO
723 : END DO
724 : END DO
725 :
726 520 : SELECT CASE (mp2_env%ri_grad%z_solver_method)
727 : CASE (z_solver_pople)
728 : CALL solve_z_vector_pople(qs_env, mp2_env, unit_nr, &
729 : mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
730 248 : L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
731 : CASE (z_solver_cg, z_solver_richardson, z_solver_sd)
732 : CALL solve_z_vector_cg(qs_env, mp2_env, unit_nr, &
733 : mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
734 24 : L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
735 : CASE DEFAULT
736 272 : CPABORT("Unknown solver")
737 : END SELECT
738 :
739 272 : CALL cp_fm_release(precond)
740 :
741 272 : CALL timestop(handle)
742 :
743 544 : END SUBROUTINE solve_z_vector_eq_low
744 :
745 : ! **************************************************************************************************
746 : !> \brief ...
747 : !> \param qs_env ...
748 : !> \param mp2_env ...
749 : !> \param unit_nr ...
750 : !> \param mo_coeff_o ...
751 : !> \param mo_coeff_v ...
752 : !> \param Eigenval ...
753 : !> \param p_env ...
754 : !> \param L_jb ...
755 : !> \param fm_G_mu_nu ...
756 : !> \param fm_back ...
757 : !> \param P_ia ...
758 : !> \param precond ...
759 : ! **************************************************************************************************
760 248 : SUBROUTINE solve_z_vector_pople(qs_env, mp2_env, unit_nr, &
761 496 : mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
762 248 : L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
763 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
764 : TYPE(mp2_type), INTENT(IN) :: mp2_env
765 : INTEGER, INTENT(IN) :: unit_nr
766 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff_o, mo_coeff_v
767 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
768 : TYPE(qs_p_env_type), INTENT(IN), POINTER :: p_env
769 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: L_jb
770 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_G_mu_nu
771 : TYPE(cp_fm_type), INTENT(IN) :: fm_back
772 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: P_ia, precond
773 :
774 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_z_vector_pople'
775 :
776 : INTEGER :: cycle_counter, handle, iiB, iiter, &
777 : ispin, max_num_iter, nspins, &
778 : transf_type_in
779 : LOGICAL :: converged
780 : REAL(KIND=dp) :: conv, eps_conv, scale_cphf, t1, t2
781 248 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: proj_bi_xj, temp_vals, x_norm, xi_b
782 248 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: A_small, b_small, xi_Axi
783 : TYPE(cp_fm_struct_p_type), ALLOCATABLE, &
784 248 : DIMENSION(:) :: fm_struct_tmp
785 248 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: b_i, residual
786 248 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: Ax, xn
787 :
788 248 : CALL timeset(routineN, handle)
789 :
790 248 : nspins = SIZE(eigenval, 2)
791 :
792 248 : eps_conv = mp2_env%ri_grad%cphf_eps_conv
793 :
794 248 : IF (accurate_dot_product_spin(L_jb, L_jb) >= (eps_conv*eps_conv)) THEN
795 :
796 248 : max_num_iter = mp2_env%ri_grad%cphf_max_num_iter
797 248 : scale_cphf = mp2_env%ri_grad%scale_step_size
798 :
799 : ! set the transformation type (equal for all methods all updates)
800 248 : transf_type_in = 3
801 :
802 : ! set convergence flag
803 248 : converged = .FALSE.
804 :
805 2490 : ALLOCATE (fm_struct_tmp(nspins), b_i(nspins), residual(nspins))
806 582 : DO ispin = 1, nspins
807 334 : fm_struct_tmp(ispin)%struct => P_ia(ispin)%matrix_struct
808 :
809 334 : CALL cp_fm_create(b_i(ispin), fm_struct_tmp(ispin)%struct, name="b_i")
810 334 : CALL cp_fm_set_all(b_i(ispin), 0.0_dp)
811 14656 : b_i(ispin)%local_data(:, :) = precond(ispin)%local_data(:, :)*L_jb(ispin)%local_data(:, :)
812 :
813 : ! create the residual vector (r), we check convergence on the norm of
814 : ! this vector r=(Ax-b)
815 334 : CALL cp_fm_create(residual(ispin), fm_struct_tmp(ispin)%struct, name="residual")
816 582 : CALL cp_fm_set_all(residual(ispin), 0.0_dp)
817 : END DO
818 :
819 248 : IF (unit_nr > 0) THEN
820 124 : WRITE (unit_nr, *)
821 124 : WRITE (unit_nr, '(T3,A)') "MP2_CPHF| Iterative solution of Z-Vector equations (Pople's method)"
822 124 : WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', eps_conv
823 124 : WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', max_num_iter
824 124 : WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Scaling of initial guess: ', scale_cphf
825 124 : WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
826 124 : WRITE (unit_nr, '(T4,A,T15,A,T33,A)') 'Step', 'Time', 'Convergence'
827 124 : WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
828 : END IF
829 :
830 11882 : ALLOCATE (xn(nspins, max_num_iter))
831 11634 : ALLOCATE (Ax(nspins, max_num_iter))
832 744 : ALLOCATE (x_norm(max_num_iter))
833 496 : ALLOCATE (xi_b(max_num_iter))
834 992 : ALLOCATE (xi_Axi(max_num_iter, 0:max_num_iter))
835 4946 : x_norm = 0.0_dp
836 4946 : xi_b = 0.0_dp
837 139110 : xi_Axi = 0.0_dp
838 :
839 248 : cycle_counter = 0
840 1056 : DO iiter = 1, max_num_iter
841 1034 : cycle_counter = cycle_counter + 1
842 :
843 1034 : t1 = m_walltime()
844 :
845 : ! create and update x_i (orthogonalization with previous vectors)
846 2446 : DO ispin = 1, nspins
847 1412 : CALL cp_fm_create(xn(ispin, iiter), fm_struct_tmp(ispin)%struct, name="xi")
848 2446 : CALL cp_fm_set_all(xn(ispin, iiter), 0.0_dp)
849 : END DO
850 :
851 2854 : ALLOCATE (proj_bi_xj(iiter - 1))
852 3028 : proj_bi_xj = 0.0_dp
853 : ! first compute the projection of the actual b_i into all previous x_i
854 : ! already scaled with the norm of each x_i
855 3028 : DO iiB = 1, iiter - 1
856 3028 : proj_bi_xj(iiB) = proj_bi_xj(iiB) + accurate_dot_product_spin(b_i, xn(:, iiB))/x_norm(iiB)
857 : END DO
858 :
859 : ! update actual x_i
860 2446 : DO ispin = 1, nspins
861 67285 : xn(ispin, iiter)%local_data(:, :) = scale_cphf*b_i(ispin)%local_data(:, :)
862 5338 : DO iiB = 1, iiter - 1
863 : xn(ispin, iiter)%local_data(:, :) = xn(ispin, iiter)%local_data(:, :) - &
864 151476 : xn(ispin, iiB)%local_data(:, :)*proj_bi_xj(iiB)
865 : END DO
866 : END DO
867 1034 : DEALLOCATE (proj_bi_xj)
868 :
869 : ! create Ax(iiter) that will store the matrix vector product for this cycle
870 2446 : DO ispin = 1, nspins
871 1412 : CALL cp_fm_create(Ax(ispin, iiter), fm_struct_tmp(ispin)%struct, name="Ai")
872 2446 : CALL cp_fm_set_all(Ax(ispin, iiter), 0.0_dp)
873 : END DO
874 :
875 : CALL cphf_like_update(qs_env, mo_coeff_o, &
876 : mo_coeff_v, Eigenval, p_env, &
877 : xn(:, iiter), fm_G_mu_nu, fm_back, transf_type_in, &
878 1034 : Ax(:, iiter))
879 :
880 : ! in order to reduce the number of parallel sums here we
881 : ! cluster all necessary scalar products into a single vector
882 : ! temp_vals contains:
883 : ! 1:iiter -> <Ax_i|x_j>
884 : ! iiter+1 -> <x_i|b>
885 : ! iiter+2 -> <x_i|x_i>
886 :
887 3102 : ALLOCATE (temp_vals(iiter + 2))
888 6130 : temp_vals = 0.0_dp
889 : ! <Ax_i|x_j>
890 4062 : DO iiB = 1, iiter
891 4062 : temp_vals(iiB) = temp_vals(iiB) + accurate_dot_product_spin(Ax(:, iiter), xn(:, iiB))
892 : END DO
893 : ! <x_i|b>
894 1034 : temp_vals(iiter + 1) = temp_vals(iiter + 1) + accurate_dot_product_spin(xn(:, iiter), L_jb)
895 : ! norm
896 1034 : temp_vals(iiter + 2) = temp_vals(iiter + 2) + accurate_dot_product_spin(xn(:, iiter), xn(:, iiter))
897 : ! update <Ax_i|x_j>, <x_i|b> and norm <x_i|x_i>
898 4062 : xi_Axi(iiter, 1:iiter) = temp_vals(1:iiter)
899 4062 : xi_Axi(1:iiter, iiter) = temp_vals(1:iiter)
900 1034 : xi_b(iiter) = temp_vals(iiter + 1)
901 1034 : x_norm(iiter) = temp_vals(iiter + 2)
902 1034 : DEALLOCATE (temp_vals)
903 :
904 : ! solve reduced system
905 1034 : IF (ALLOCATED(A_small)) DEALLOCATE (A_small)
906 1034 : IF (ALLOCATED(b_small)) DEALLOCATE (b_small)
907 4136 : ALLOCATE (A_small(iiter, iiter))
908 3102 : ALLOCATE (b_small(iiter, 1))
909 16110 : A_small(1:iiter, 1:iiter) = xi_Axi(1:iiter, 1:iiter)
910 4062 : b_small(1:iiter, 1) = xi_b(1:iiter)
911 :
912 1034 : CALL solve_system(matrix=A_small, mysize=iiter, eigenvectors=b_small)
913 :
914 : ! check for convergence
915 2446 : DO ispin = 1, nspins
916 1412 : CALL cp_fm_set_all(residual(ispin), 0.0_dp)
917 5716 : DO iiB = 1, iiter
918 : residual(ispin)%local_data(:, :) = &
919 : residual(ispin)%local_data(:, :) + &
920 218761 : b_small(iiB, 1)*Ax(ispin, iiB)%local_data(:, :)
921 : END DO
922 :
923 : residual(ispin)%local_data(:, :) = &
924 : residual(ispin)%local_data(:, :) - &
925 68319 : L_jb(ispin)%local_data(:, :)
926 : END DO
927 :
928 1034 : conv = SQRT(accurate_dot_product_spin(residual, residual))
929 :
930 1034 : t2 = m_walltime()
931 :
932 1034 : IF (unit_nr > 0) THEN
933 517 : WRITE (unit_nr, '(T3,I5,T13,F6.1,11X,F14.8)') iiter, t2 - t1, conv
934 517 : CALL m_flush(unit_nr)
935 : END IF
936 :
937 1034 : IF (conv <= eps_conv) THEN
938 : converged = .TRUE.
939 : EXIT
940 : END IF
941 :
942 : ! update b_i for the next round
943 1926 : DO ispin = 1, nspins
944 : b_i(ispin)%local_data(:, :) = b_i(ispin)%local_data(:, :) &
945 : + precond(ispin)%local_data(:, :) &
946 54956 : *Ax(ispin, iiter)%local_data(:, :)
947 : END DO
948 :
949 830 : scale_cphf = 1.0_dp
950 :
951 : END DO
952 :
953 248 : IF (unit_nr > 0) THEN
954 124 : WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
955 124 : IF (converged) THEN
956 113 : WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations converged in', cycle_counter, ' steps'
957 : ELSE
958 11 : WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations NOT converged in', cycle_counter, ' steps'
959 : END IF
960 : END IF
961 :
962 : ! store solution into P_ia
963 1282 : DO iiter = 1, cycle_counter
964 2694 : DO ispin = 1, nspins
965 : P_ia(ispin)%local_data(:, :) = P_ia(ispin)%local_data(:, :) + &
966 68319 : b_small(iiter, 1)*xn(ispin, iiter)%local_data(:, :)
967 : END DO
968 : END DO
969 :
970 : ! Release arrays
971 248 : DEALLOCATE (x_norm)
972 248 : DEALLOCATE (xi_b)
973 248 : DEALLOCATE (xi_Axi)
974 :
975 248 : CALL cp_fm_release(b_i)
976 248 : CALL cp_fm_release(residual)
977 248 : CALL cp_fm_release(Ax)
978 248 : CALL cp_fm_release(xn)
979 :
980 : ELSE
981 0 : IF (unit_nr > 0) THEN
982 0 : WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
983 0 : WRITE (unit_nr, '(T3,A)') 'Residual smaller than EPS_CONV. Skip solution of Z-vector equation.'
984 : END IF
985 : END IF
986 :
987 248 : CALL timestop(handle)
988 :
989 248 : END SUBROUTINE solve_z_vector_pople
990 :
991 : ! **************************************************************************************************
992 : !> \brief ...
993 : !> \param qs_env ...
994 : !> \param mp2_env ...
995 : !> \param unit_nr ...
996 : !> \param mo_coeff_o ...
997 : !> \param mo_coeff_v ...
998 : !> \param Eigenval ...
999 : !> \param p_env ...
1000 : !> \param L_jb ...
1001 : !> \param fm_G_mu_nu ...
1002 : !> \param fm_back ...
1003 : !> \param P_ia ...
1004 : !> \param precond ...
1005 : ! **************************************************************************************************
1006 24 : SUBROUTINE solve_z_vector_cg(qs_env, mp2_env, unit_nr, &
1007 48 : mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
1008 48 : L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
1009 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1010 : TYPE(mp2_type), INTENT(IN) :: mp2_env
1011 : INTEGER, INTENT(IN) :: unit_nr
1012 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_coeff_o, mo_coeff_v
1013 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: Eigenval
1014 : TYPE(qs_p_env_type), INTENT(IN), POINTER :: p_env
1015 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: L_jb
1016 : TYPE(cp_fm_type), INTENT(INOUT) :: fm_G_mu_nu
1017 : TYPE(cp_fm_type), INTENT(IN) :: fm_back
1018 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: P_ia, precond
1019 :
1020 : CHARACTER(LEN=*), PARAMETER :: routineN = 'solve_z_vector_cg'
1021 :
1022 : INTEGER :: cycles_passed, handle, iiter, ispin, max_num_iter, nspins, restart_counter, &
1023 : restart_every, transf_type_in, z_solver_method
1024 : LOGICAL :: converged, do_restart
1025 : REAL(KIND=dp) :: eps_conv, norm_residual, norm_residual_old, &
1026 : residual_dot_diff_search_vec_new, residual_dot_diff_search_vec_old, &
1027 : residual_dot_search_vec, residual_new_dot_diff_search_vec_old, scale_result, &
1028 : scale_search, scale_step_size, search_vec_dot_A_search_vec, t1, t2
1029 : TYPE(cp_fm_struct_p_type), ALLOCATABLE, &
1030 24 : DIMENSION(:) :: fm_struct_tmp
1031 24 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: A_dot_search_vector, diff_search_vector, &
1032 24 : residual, search_vector
1033 :
1034 24 : CALL timeset(routineN, handle)
1035 :
1036 24 : max_num_iter = mp2_env%ri_grad%cphf_max_num_iter
1037 24 : eps_conv = mp2_env%ri_grad%cphf_eps_conv
1038 24 : z_solver_method = mp2_env%ri_grad%z_solver_method
1039 24 : restart_every = mp2_env%ri_grad%cphf_restart
1040 24 : scale_step_size = mp2_env%ri_grad%scale_step_size
1041 24 : transf_type_in = 3
1042 24 : nspins = SIZE(eigenval, 2)
1043 :
1044 24 : IF (unit_nr > 0) THEN
1045 12 : WRITE (unit_nr, *)
1046 8 : SELECT CASE (z_solver_method)
1047 : CASE (z_solver_cg)
1048 8 : IF (mp2_env%ri_grad%polak_ribiere) THEN
1049 2 : WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (CG with Polak-Ribiere step)'
1050 : ELSE
1051 6 : WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (CG with Fletcher-Reeves step)'
1052 : END IF
1053 : CASE (z_solver_richardson)
1054 2 : WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (Richardson method)'
1055 : CASE (z_solver_sd)
1056 2 : WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (Steepest Descent method)'
1057 : CASE DEFAULT
1058 12 : CPABORT("Unknown solver")
1059 : END SELECT
1060 12 : WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', eps_conv
1061 12 : WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', max_num_iter
1062 12 : WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Number of steps for restart: ', restart_every
1063 12 : WRITE (unit_nr, '(T3, A)') 'MP2_CPHF| Restart after no decrease'
1064 12 : WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Scaling factor of each step: ', scale_step_size
1065 12 : WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
1066 12 : WRITE (unit_nr, '(T4,A,T13,A,T28,A,T43,A)') 'Step', 'Restart', 'Time', 'Convergence'
1067 12 : WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
1068 : END IF
1069 :
1070 0 : ALLOCATE (fm_struct_tmp(nspins), residual(nspins), diff_search_vector(nspins), &
1071 312 : search_vector(nspins), A_dot_search_vector(nspins))
1072 48 : DO ispin = 1, nspins
1073 24 : fm_struct_tmp(ispin)%struct => P_ia(ispin)%matrix_struct
1074 :
1075 24 : CALL cp_fm_create(residual(ispin), fm_struct_tmp(ispin)%struct, name="residual")
1076 24 : CALL cp_fm_set_all(residual(ispin), 0.0_dp)
1077 :
1078 24 : CALL cp_fm_create(diff_search_vector(ispin), fm_struct_tmp(ispin)%struct, name="difference search vector")
1079 24 : CALL cp_fm_set_all(diff_search_vector(ispin), 0.0_dp)
1080 :
1081 24 : CALL cp_fm_create(search_vector(ispin), fm_struct_tmp(ispin)%struct, name="search vector")
1082 24 : CALL cp_fm_set_all(search_vector(ispin), 0.0_dp)
1083 :
1084 24 : CALL cp_fm_create(A_dot_search_vector(ispin), fm_struct_tmp(ispin)%struct, name="A times search vector")
1085 48 : CALL cp_fm_set_all(A_dot_search_vector(ispin), 0.0_dp)
1086 : END DO
1087 :
1088 24 : converged = .FALSE.
1089 24 : cycles_passed = max_num_iter
1090 : ! By that, we enforce the setup of the matrices
1091 24 : do_restart = .TRUE.
1092 :
1093 24 : t1 = m_walltime()
1094 :
1095 164 : DO iiter = 1, max_num_iter
1096 :
1097 : ! During the first iteration, P_ia=0 such that the application of the 2nd order matrix is zero
1098 160 : IF (do_restart) THEN
1099 : ! We do not consider the first step to be a restart
1100 : ! Do not recalculate residual if it is already enforced to save FLOPs
1101 52 : IF (.NOT. mp2_env%ri_grad%recalc_residual .OR. (iiter == 1)) THEN
1102 52 : IF (iiter > 1) THEN
1103 : CALL cphf_like_update(qs_env, mo_coeff_o, &
1104 : mo_coeff_v, Eigenval, p_env, &
1105 : P_ia, fm_G_mu_nu, fm_back, transf_type_in, &
1106 28 : residual)
1107 : ELSE
1108 24 : do_restart = .FALSE.
1109 :
1110 48 : DO ispin = 1, nspins
1111 48 : CALL cp_fm_set_all(residual(ispin), 0.0_dp)
1112 : END DO
1113 : END IF
1114 :
1115 104 : DO ispin = 1, nspins
1116 : residual(ispin)%local_data(:, :) = L_jb(ispin)%local_data(:, :) &
1117 3068 : - residual(ispin)%local_data(:, :)
1118 : END DO
1119 : END IF
1120 :
1121 104 : DO ispin = 1, nspins
1122 : diff_search_vector(ispin)%local_data(:, :) = &
1123 3016 : precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1124 3068 : search_vector(ispin)%local_data(:, :) = diff_search_vector(ispin)%local_data(:, :)
1125 : END DO
1126 :
1127 : restart_counter = 1
1128 : END IF
1129 :
1130 160 : norm_residual_old = SQRT(accurate_dot_product_spin(residual, residual))
1131 :
1132 160 : residual_dot_diff_search_vec_old = accurate_dot_product_spin(residual, diff_search_vector)
1133 :
1134 : CALL cphf_like_update(qs_env, mo_coeff_o, &
1135 : mo_coeff_v, Eigenval, p_env, &
1136 : search_vector, fm_G_mu_nu, fm_back, transf_type_in, &
1137 160 : A_dot_search_vector)
1138 :
1139 160 : IF (z_solver_method /= z_solver_richardson) THEN
1140 120 : search_vec_dot_A_search_vec = accurate_dot_product_spin(search_vector, A_dot_search_vector)
1141 :
1142 120 : IF (z_solver_method == z_solver_cg) THEN
1143 94 : scale_result = residual_dot_diff_search_vec_old/search_vec_dot_A_search_vec
1144 : ELSE
1145 26 : residual_dot_search_vec = accurate_dot_product_spin(residual, search_vector)
1146 26 : scale_result = residual_dot_search_vec/search_vec_dot_A_search_vec
1147 : END IF
1148 :
1149 120 : scale_result = scale_result*scale_step_size
1150 :
1151 : ELSE
1152 :
1153 : scale_result = scale_step_size
1154 :
1155 : END IF
1156 :
1157 320 : DO ispin = 1, nspins
1158 : P_ia(ispin)%local_data(:, :) = P_ia(ispin)%local_data(:, :) &
1159 9440 : + scale_result*search_vector(ispin)%local_data(:, :)
1160 : END DO
1161 :
1162 160 : IF (.NOT. mp2_env%ri_grad%recalc_residual) THEN
1163 :
1164 284 : DO ispin = 1, nspins
1165 : residual(ispin)%local_data(:, :) = residual(ispin)%local_data(:, :) &
1166 8378 : - scale_result*A_dot_search_vector(ispin)%local_data(:, :)
1167 : END DO
1168 : ELSE
1169 : CALL cphf_like_update(qs_env, mo_coeff_o, &
1170 : mo_coeff_v, Eigenval, p_env, &
1171 : P_ia, fm_G_mu_nu, fm_back, transf_type_in, &
1172 18 : residual)
1173 :
1174 36 : DO ispin = 1, nspins
1175 1062 : residual(ispin)%local_data(:, :) = L_jb(ispin)%local_data(:, :) - residual(ispin)%local_data(:, :)
1176 : END DO
1177 : END IF
1178 :
1179 160 : norm_residual = SQRT(accurate_dot_product_spin(residual, residual))
1180 :
1181 160 : t2 = m_walltime()
1182 :
1183 160 : IF (unit_nr > 0) THEN
1184 80 : WRITE (unit_nr, '(T3,I4,T16,L1,T26,F6.1,8X,F14.8)') iiter, do_restart, t2 - t1, norm_residual
1185 80 : CALL m_flush(unit_nr)
1186 : END IF
1187 :
1188 160 : IF (norm_residual <= eps_conv) THEN
1189 20 : converged = .TRUE.
1190 20 : cycles_passed = iiter
1191 20 : EXIT
1192 : END IF
1193 :
1194 140 : t1 = m_walltime()
1195 :
1196 140 : IF (z_solver_method == z_solver_richardson) THEN
1197 80 : DO ispin = 1, nspins
1198 : search_vector(ispin)%local_data(:, :) = &
1199 2360 : scale_step_size*precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1200 : END DO
1201 100 : ELSE IF (z_solver_method == z_solver_sd) THEN
1202 44 : DO ispin = 1, nspins
1203 : search_vector(ispin)%local_data(:, :) = &
1204 1298 : precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1205 : END DO
1206 : ELSE
1207 78 : IF (mp2_env%ri_grad%polak_ribiere) &
1208 28 : residual_new_dot_diff_search_vec_old = accurate_dot_product_spin(residual, diff_search_vector)
1209 :
1210 156 : DO ispin = 1, nspins
1211 : diff_search_vector(ispin)%local_data(:, :) = &
1212 4602 : precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1213 : END DO
1214 :
1215 78 : residual_dot_diff_search_vec_new = accurate_dot_product_spin(residual, diff_search_vector)
1216 :
1217 78 : scale_search = residual_dot_diff_search_vec_new/residual_dot_diff_search_vec_old
1218 78 : IF (mp2_env%ri_grad%polak_ribiere) scale_search = &
1219 28 : scale_search - residual_new_dot_diff_search_vec_old/residual_dot_diff_search_vec_old
1220 :
1221 156 : DO ispin = 1, nspins
1222 : search_vector(ispin)%local_data(:, :) = scale_search*search_vector(ispin)%local_data(:, :) &
1223 4602 : + diff_search_vector(ispin)%local_data(:, :)
1224 : END DO
1225 :
1226 : ! Make new to old
1227 140 : residual_dot_diff_search_vec_old = residual_dot_diff_search_vec_new
1228 : END IF
1229 :
1230 : ! Check whether the residual decrease or restart is enforced and ask for restart
1231 140 : do_restart = (norm_residual >= norm_residual_old .OR. (MOD(restart_counter, restart_every) == 0))
1232 :
1233 140 : restart_counter = restart_counter + 1
1234 144 : norm_residual_old = norm_residual
1235 :
1236 : END DO
1237 :
1238 24 : IF (unit_nr > 0) THEN
1239 12 : WRITE (unit_nr, '(T4,A)') REPEAT("-", 40)
1240 12 : IF (converged) THEN
1241 10 : WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations converged in', cycles_passed, ' steps'
1242 : ELSE
1243 2 : WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations NOT converged in', max_num_iter, ' steps'
1244 : END IF
1245 : END IF
1246 :
1247 24 : DEALLOCATE (fm_struct_tmp)
1248 24 : CALL cp_fm_release(residual)
1249 24 : CALL cp_fm_release(diff_search_vector)
1250 24 : CALL cp_fm_release(search_vector)
1251 24 : CALL cp_fm_release(A_dot_search_vector)
1252 :
1253 24 : CALL timestop(handle)
1254 :
1255 48 : END SUBROUTINE solve_z_vector_cg
1256 :
1257 : ! **************************************************************************************************
1258 : !> \brief ...
1259 : !> \param matrix1 ...
1260 : !> \param matrix2 ...
1261 : !> \return ...
1262 : ! **************************************************************************************************
1263 9104 : FUNCTION accurate_dot_product_spin(matrix1, matrix2) RESULT(dotproduct)
1264 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: matrix1, matrix2
1265 : REAL(KIND=dp) :: dotproduct
1266 :
1267 : INTEGER :: ispin
1268 :
1269 9104 : dotproduct = 0.0_dp
1270 21602 : DO ispin = 1, SIZE(matrix1)
1271 21602 : dotproduct = dotproduct + accurate_dot_product(matrix1(ispin)%local_data, matrix2(ispin)%local_data)
1272 : END DO
1273 9104 : CALL matrix1(1)%matrix_struct%para_env%sum(dotproduct)
1274 :
1275 9104 : END FUNCTION accurate_dot_product_spin
1276 :
1277 : ! **************************************************************************************************
1278 : !> \brief ...
1279 : !> \param qs_env ...
1280 : ! **************************************************************************************************
1281 272 : SUBROUTINE update_mp2_forces(qs_env)
1282 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1283 :
1284 : CHARACTER(LEN=*), PARAMETER :: routineN = 'update_mp2_forces'
1285 :
1286 : INTEGER :: alpha, beta, handle, idir, iounit, &
1287 : ispin, nimages, nocc, nspins
1288 : INTEGER, DIMENSION(3) :: comp
1289 : LOGICAL :: do_exx, do_hfx, use_virial
1290 : REAL(KIND=dp) :: e_dummy, e_hartree, e_xc, ehartree, exc
1291 : REAL(KIND=dp), DIMENSION(3) :: deb
1292 : REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_virial
1293 : TYPE(admm_type), POINTER :: admm_env
1294 : TYPE(cell_type), POINTER :: cell
1295 : TYPE(cp_logger_type), POINTER :: logger
1296 : TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:), &
1297 272 : TARGET :: matrix_ks_aux
1298 272 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p_mp2, &
1299 272 : matrix_p_mp2_admm, matrix_s, rho1, &
1300 272 : rho_ao, rho_ao_aux, scrm
1301 272 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p, rho_ao_kp, scrm_kp
1302 : TYPE(dft_control_type), POINTER :: dft_control
1303 272 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1304 : TYPE(linres_control_type), POINTER :: linres_control
1305 272 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1306 : TYPE(mp_para_env_type), POINTER :: para_env
1307 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1308 272 : POINTER :: sab_orb
1309 : TYPE(pw_c1d_gs_type) :: pot_g, rho_tot_g, temp_pw_g
1310 272 : TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: dvg
1311 272 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_mp2_g
1312 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
1313 : TYPE(pw_env_type), POINTER :: pw_env
1314 : TYPE(pw_poisson_type), POINTER :: poisson_env
1315 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1316 : TYPE(pw_r3d_rs_type) :: pot_r, vh_rspace, vhxc_rspace
1317 272 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_mp2_r, rho_mp2_r_aux, rho_r, &
1318 272 : tau_mp2_r, vadmm_rspace, vtau_rspace, &
1319 272 : vxc_rspace
1320 : TYPE(qs_dispersion_type), POINTER :: dispersion_env
1321 : TYPE(qs_energy_type), POINTER :: energy
1322 272 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1323 : TYPE(qs_ks_env_type), POINTER :: ks_env
1324 : TYPE(qs_p_env_type), POINTER :: p_env
1325 : TYPE(qs_rho_type), POINTER :: rho, rho_aux
1326 : TYPE(section_vals_type), POINTER :: hfx_section, hfx_sections, input, &
1327 : xc_section
1328 : TYPE(task_list_type), POINTER :: task_list_aux_fit
1329 : TYPE(virial_type), POINTER :: virial
1330 :
1331 272 : CALL timeset(routineN, handle)
1332 :
1333 272 : NULLIFY (input, pw_env, matrix_s, rho, energy, force, virial, &
1334 272 : matrix_p_mp2, matrix_p_mp2_admm, matrix_ks, rho_core)
1335 : CALL get_qs_env(qs_env, &
1336 : ks_env=ks_env, &
1337 : dft_control=dft_control, &
1338 : pw_env=pw_env, &
1339 : input=input, &
1340 : mos=mos, &
1341 : para_env=para_env, &
1342 : matrix_s=matrix_s, &
1343 : matrix_ks=matrix_ks, &
1344 : matrix_p_mp2=matrix_p_mp2, &
1345 : matrix_p_mp2_admm=matrix_p_mp2_admm, &
1346 : rho=rho, &
1347 : cell=cell, &
1348 : force=force, &
1349 : virial=virial, &
1350 : sab_orb=sab_orb, &
1351 : energy=energy, &
1352 : rho_core=rho_core, &
1353 272 : x_data=x_data)
1354 :
1355 272 : logger => cp_get_default_logger()
1356 : iounit = cp_print_key_unit_nr(logger, input, "DFT%XC%WF_CORRELATION%PRINT", &
1357 272 : extension=".mp2Log")
1358 :
1359 272 : do_exx = .FALSE.
1360 272 : IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
1361 24 : hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
1362 24 : CALL section_vals_get(hfx_section, explicit=do_exx)
1363 : END IF
1364 :
1365 272 : nimages = dft_control%nimages
1366 272 : CPASSERT(nimages == 1)
1367 :
1368 272 : p_env => qs_env%mp2_env%ri_grad%p_env
1369 :
1370 272 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp, rho_r=rho_r, rho_g=rho_g)
1371 272 : nspins = SIZE(rho_ao)
1372 :
1373 : ! check if we have to calculate the virial
1374 272 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1375 272 : IF (use_virial) virial%pv_calculate = .TRUE.
1376 :
1377 272 : CALL zero_qs_force(force)
1378 272 : IF (use_virial) CALL zero_virial(virial, .FALSE.)
1379 :
1380 630 : DO ispin = 1, nspins
1381 630 : CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
1382 : END DO
1383 :
1384 : ! kinetic energy
1385 272 : NULLIFY (scrm_kp)
1386 272 : matrix_p(1:nspins, 1:1) => rho_ao(1:nspins)
1387 : CALL kinetic_energy_matrix(qs_env, matrixkp_t=scrm_kp, matrix_p=matrix_p, &
1388 : calculate_forces=.TRUE., &
1389 272 : debug_forces=debug_forces, debug_stress=debug_forces)
1390 272 : CALL dbcsr_deallocate_matrix_set(scrm_kp)
1391 :
1392 272 : scrm_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
1393 : CALL core_matrices(qs_env, scrm_kp, rho_ao_kp, .TRUE., 1, &
1394 272 : debug_forces=debug_forces, debug_stress=debug_forces)
1395 :
1396 : ! Get the different components of the KS potential
1397 272 : NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)
1398 272 : IF (use_virial) THEN
1399 50 : h_stress = 0.0_dp
1400 50 : CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
1401 : ! Update virial
1402 650 : virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
1403 650 : virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
1404 50 : IF (.NOT. do_exx) THEN
1405 650 : virial%pv_exc = virial%pv_exc - virial%pv_xc
1406 650 : virial%pv_virial = virial%pv_virial - virial%pv_xc
1407 : ELSE
1408 0 : virial%pv_xc = 0.0_dp
1409 : END IF
1410 : IF (debug_forces) THEN
1411 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: xc ", third_tr(h_stress)
1412 50 : CALL para_env%sum(virial%pv_xc(1, 1))
1413 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Corr xc ", third_tr(virial%pv_xc)
1414 : END IF
1415 : ELSE
1416 222 : CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc)
1417 : END IF
1418 :
1419 : ! Vhxc
1420 272 : CALL get_qs_env(qs_env, pw_env=pw_env)
1421 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1422 272 : poisson_env=poisson_env)
1423 272 : CALL auxbas_pw_pool%create_pw(vhxc_rspace)
1424 872 : IF (use_virial) h_stress = virial%pv_virial
1425 : IF (debug_forces) THEN
1426 1088 : deb(1:3) = force(1)%rho_elec(1:3, 1)
1427 272 : IF (use_virial) e_dummy = third_tr(h_stress)
1428 : END IF
1429 272 : IF (do_exx) THEN
1430 16 : DO ispin = 1, nspins
1431 8 : CALL pw_transfer(vh_rspace, vhxc_rspace)
1432 8 : CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
1433 : CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1434 : hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
1435 8 : qs_env=qs_env, calculate_forces=.TRUE.)
1436 8 : CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
1437 8 : CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1438 : CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1439 : hmat=matrix_ks(ispin), pmat=matrix_p_mp2(ispin), &
1440 8 : qs_env=qs_env, calculate_forces=.TRUE.)
1441 16 : IF (ASSOCIATED(vtau_rspace)) THEN
1442 : CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
1443 : hmat=matrix_ks(ispin), pmat=matrix_p_mp2(ispin), &
1444 0 : qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.TRUE.)
1445 : END IF
1446 : END DO
1447 : ELSE
1448 614 : DO ispin = 1, nspins
1449 350 : CALL pw_transfer(vh_rspace, vhxc_rspace)
1450 350 : CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1451 : CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1452 : hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
1453 350 : qs_env=qs_env, calculate_forces=.TRUE.)
1454 614 : IF (ASSOCIATED(vtau_rspace)) THEN
1455 : CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
1456 : hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
1457 60 : qs_env=qs_env, calculate_forces=.TRUE., compute_tau=.TRUE.)
1458 : END IF
1459 : END DO
1460 : END IF
1461 : IF (debug_forces) THEN
1462 1088 : deb(1:3) = force(1)%rho_elec(1:3, 1) - deb(1:3)
1463 272 : CALL para_env%sum(deb)
1464 272 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dVhxc ", deb
1465 272 : IF (use_virial) THEN
1466 50 : e_dummy = third_tr(virial%pv_virial) - e_dummy
1467 50 : CALL para_env%sum(e_dummy)
1468 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Vhxc ", e_dummy
1469 : END IF
1470 : END IF
1471 272 : IF (use_virial) THEN
1472 650 : h_stress = virial%pv_virial - h_stress
1473 650 : virial%pv_ehartree = virial%pv_ehartree + h_stress
1474 :
1475 50 : CALL qs_rho_get(p_env%rho1, rho_r=rho_mp2_r, tau_r=tau_mp2_r)
1476 50 : e_xc = 0.0_dp
1477 124 : DO ispin = 1, nspins
1478 : ! The potentials have been scaled in ks_ref_potential, but for pw_integral_ab, we need the unscaled potentials
1479 74 : CALL pw_scale(vxc_rspace(ispin), 1.0_dp/vxc_rspace(ispin)%pw_grid%dvol)
1480 74 : e_xc = e_xc + pw_integral_ab(rho_mp2_r(ispin), vxc_rspace(ispin))
1481 74 : IF (ASSOCIATED(vtau_rspace)) CALL pw_scale(vtau_rspace(ispin), 1.0_dp/vtau_rspace(ispin)%pw_grid%dvol)
1482 124 : IF (ASSOCIATED(vtau_rspace)) e_xc = e_xc + pw_integral_ab(tau_mp2_r(ispin), vtau_rspace(ispin))
1483 : END DO
1484 50 : IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG VIRIAL:: vxc*d1 ", e_xc
1485 200 : DO alpha = 1, 3
1486 150 : virial%pv_exc(alpha, alpha) = virial%pv_exc(alpha, alpha) - e_xc/REAL(para_env%num_pe, dp)
1487 200 : virial%pv_virial(alpha, alpha) = virial%pv_virial(alpha, alpha) - e_xc/REAL(para_env%num_pe, dp)
1488 : END DO
1489 : END IF
1490 630 : DO ispin = 1, nspins
1491 358 : CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
1492 630 : IF (ASSOCIATED(vtau_rspace)) THEN
1493 60 : CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
1494 : END IF
1495 : END DO
1496 272 : DEALLOCATE (vxc_rspace)
1497 272 : CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
1498 272 : IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
1499 :
1500 630 : DO ispin = 1, nspins
1501 630 : CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
1502 : END DO
1503 :
1504 : ! pw stuff
1505 272 : NULLIFY (poisson_env, auxbas_pw_pool)
1506 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1507 272 : poisson_env=poisson_env)
1508 :
1509 : ! get some of the grids ready
1510 272 : CALL auxbas_pw_pool%create_pw(pot_r)
1511 272 : CALL auxbas_pw_pool%create_pw(pot_g)
1512 272 : CALL auxbas_pw_pool%create_pw(rho_tot_g)
1513 :
1514 272 : CALL pw_zero(rho_tot_g)
1515 :
1516 272 : CALL qs_rho_get(p_env%rho1, rho_r=rho_mp2_r, rho_g=rho_mp2_g)
1517 630 : DO ispin = 1, nspins
1518 630 : CALL pw_axpy(rho_mp2_g(ispin), rho_tot_g)
1519 : END DO
1520 :
1521 272 : IF (use_virial) THEN
1522 200 : ALLOCATE (dvg(3))
1523 200 : DO idir = 1, 3
1524 200 : CALL auxbas_pw_pool%create_pw(dvg(idir))
1525 : END DO
1526 50 : CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g, dvhartree=dvg)
1527 : ELSE
1528 222 : CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g)
1529 : END IF
1530 :
1531 272 : CALL pw_transfer(pot_g, pot_r)
1532 272 : CALL pw_scale(pot_r, pot_r%pw_grid%dvol)
1533 272 : CALL pw_axpy(pot_r, vh_rspace)
1534 :
1535 : ! calculate core forces
1536 272 : CALL integrate_v_core_rspace(vh_rspace, qs_env)
1537 : IF (debug_forces) THEN
1538 1088 : deb(:) = force(1)%rho_core(:, 1)
1539 272 : CALL para_env%sum(deb)
1540 272 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: core ", deb
1541 272 : IF (use_virial) THEN
1542 50 : e_dummy = third_tr(virial%pv_virial) - e_dummy
1543 50 : CALL para_env%sum(e_dummy)
1544 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: core ", e_dummy
1545 : END IF
1546 : END IF
1547 272 : CALL auxbas_pw_pool%give_back_pw(vh_rspace)
1548 :
1549 272 : IF (use_virial) THEN
1550 : ! update virial if necessary with the volume term
1551 : ! first create pw auxiliary stuff
1552 50 : CALL auxbas_pw_pool%create_pw(temp_pw_g)
1553 :
1554 : ! make a copy of the MP2 density in G space
1555 50 : CALL pw_copy(rho_tot_g, temp_pw_g)
1556 :
1557 : ! calculate total SCF density and potential
1558 50 : CALL pw_copy(rho_g(1), rho_tot_g)
1559 50 : IF (nspins == 2) CALL pw_axpy(rho_g(2), rho_tot_g)
1560 50 : CALL pw_axpy(rho_core, rho_tot_g)
1561 50 : CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g)
1562 :
1563 : ! finally update virial with the volume contribution
1564 50 : e_hartree = pw_integral_ab(temp_pw_g, pot_g)
1565 50 : IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: vh1*d0 ", e_hartree
1566 :
1567 50 : h_stress = 0.0_dp
1568 200 : DO alpha = 1, 3
1569 150 : comp = 0
1570 150 : comp(alpha) = 1
1571 150 : CALL pw_copy(pot_g, rho_tot_g)
1572 150 : CALL pw_derive(rho_tot_g, comp)
1573 150 : h_stress(alpha, alpha) = -e_hartree
1574 500 : DO beta = alpha, 3
1575 : h_stress(alpha, beta) = h_stress(alpha, beta) &
1576 300 : - 2.0_dp*pw_integral_ab(rho_tot_g, dvg(beta))/fourpi
1577 450 : h_stress(beta, alpha) = h_stress(alpha, beta)
1578 : END DO
1579 : END DO
1580 50 : IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hartree ", third_tr(h_stress)
1581 :
1582 : ! free stuff
1583 50 : CALL auxbas_pw_pool%give_back_pw(temp_pw_g)
1584 200 : DO idir = 1, 3
1585 200 : CALL auxbas_pw_pool%give_back_pw(dvg(idir))
1586 : END DO
1587 50 : DEALLOCATE (dvg)
1588 :
1589 650 : virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
1590 650 : virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
1591 :
1592 : END IF
1593 :
1594 630 : DO ispin = 1, nspins
1595 358 : CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
1596 630 : IF (dft_control%do_admm) CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
1597 : END DO
1598 :
1599 272 : CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
1600 :
1601 272 : IF (dft_control%do_admm) THEN
1602 40 : CALL get_qs_env(qs_env, admm_env=admm_env)
1603 40 : xc_section => admm_env%xc_section_primary
1604 : ELSE
1605 232 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1606 : END IF
1607 :
1608 272 : IF (use_virial) THEN
1609 50 : h_stress = 0.0_dp
1610 650 : pv_virial = virial%pv_virial
1611 : END IF
1612 : IF (debug_forces) THEN
1613 1088 : deb = force(1)%rho_elec(1:3, 1)
1614 272 : IF (use_virial) e_dummy = third_tr(pv_virial)
1615 : END IF
1616 272 : CALL apply_2nd_order_kernel(qs_env, p_env, .FALSE., .TRUE., use_virial, h_stress)
1617 272 : IF (use_virial) THEN
1618 650 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_virial)
1619 : IF (debug_forces) THEN
1620 650 : e_dummy = third_tr(virial%pv_virial - pv_virial)
1621 50 : CALL para_env%sum(e_dummy)
1622 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Kh ", e_dummy
1623 : END IF
1624 650 : virial%pv_exc = virial%pv_exc + h_stress
1625 650 : virial%pv_virial = virial%pv_virial + h_stress
1626 : IF (debug_forces) THEN
1627 50 : e_dummy = third_tr(h_stress)
1628 50 : CALL para_env%sum(e_dummy)
1629 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Kxc ", e_dummy
1630 : END IF
1631 : END IF
1632 : IF (debug_forces) THEN
1633 1088 : deb(:) = force(1)%rho_elec(:, 1) - deb
1634 272 : CALL para_env%sum(deb)
1635 272 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P0*Khxc ", deb
1636 272 : IF (use_virial) THEN
1637 50 : e_dummy = third_tr(virial%pv_virial) - e_dummy
1638 50 : CALL para_env%sum(e_dummy)
1639 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Khxc ", e_dummy
1640 : END IF
1641 : END IF
1642 :
1643 : ! hfx section
1644 272 : NULLIFY (hfx_sections)
1645 272 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1646 272 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
1647 272 : IF (do_hfx) THEN
1648 192 : IF (do_exx) THEN
1649 0 : IF (dft_control%do_admm) THEN
1650 0 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
1651 0 : CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux, rho_ao_kp=rho_ao_kp)
1652 0 : rho1 => p_env%p1_admm
1653 : ELSE
1654 0 : rho1 => p_env%p1
1655 : END IF
1656 : ELSE
1657 192 : IF (dft_control%do_admm) THEN
1658 36 : CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
1659 36 : CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux, rho_ao_kp=rho_ao_kp)
1660 90 : DO ispin = 1, nspins
1661 90 : CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
1662 : END DO
1663 36 : rho1 => p_env%p1_admm
1664 : ELSE
1665 344 : DO ispin = 1, nspins
1666 344 : CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
1667 : END DO
1668 156 : rho1 => p_env%p1
1669 : END IF
1670 : END IF
1671 :
1672 192 : IF (x_data(1, 1)%do_hfx_ri) THEN
1673 :
1674 : CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
1675 : x_data(1, 1)%general_parameter%fraction, &
1676 : rho_ao=rho_ao_kp, rho_ao_resp=rho1, &
1677 6 : use_virial=use_virial, resp_only=do_exx)
1678 :
1679 : ELSE
1680 : CALL derivatives_four_center(qs_env, rho_ao_kp, rho1, hfx_sections, para_env, &
1681 186 : 1, use_virial, resp_only=do_exx)
1682 : END IF
1683 :
1684 192 : IF (use_virial) THEN
1685 338 : virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1686 338 : virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1687 : END IF
1688 : IF (debug_forces) THEN
1689 768 : deb(1:3) = force(1)%fock_4c(1:3, 1) - deb(1:3)
1690 192 : CALL para_env%sum(deb)
1691 192 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx ", deb
1692 192 : IF (use_virial) THEN
1693 26 : e_dummy = third_tr(virial%pv_fock_4c)
1694 26 : CALL para_env%sum(e_dummy)
1695 26 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: hfx ", e_dummy
1696 : END IF
1697 : END IF
1698 :
1699 192 : IF (.NOT. do_exx) THEN
1700 192 : IF (dft_control%do_admm) THEN
1701 36 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1702 90 : DO ispin = 1, nspins
1703 90 : CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
1704 : END DO
1705 : ELSE
1706 344 : DO ispin = 1, nspins
1707 344 : CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, -1.0_dp)
1708 : END DO
1709 : END IF
1710 : END IF
1711 :
1712 192 : IF (dft_control%do_admm) THEN
1713 : IF (debug_forces) THEN
1714 144 : deb = force(1)%overlap_admm(1:3, 1)
1715 36 : IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1716 : END IF
1717 : ! The 2nd order kernel contains a factor of two in apply_xc_admm_ao which we don't need for the projection derivatives
1718 36 : IF (nspins == 1) CALL dbcsr_scale(p_env%kpp1_admm(1)%matrix, 0.5_dp)
1719 36 : CALL admm_projection_derivative(qs_env, p_env%kpp1_admm, rho_ao)
1720 : IF (debug_forces) THEN
1721 144 : deb(:) = force(1)%overlap_admm(:, 1) - deb
1722 36 : CALL para_env%sum(deb)
1723 36 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*KADMM*dS'", deb
1724 36 : IF (use_virial) THEN
1725 12 : e_dummy = third_tr(virial%pv_virial) - e_dummy
1726 12 : CALL para_env%sum(e_dummy)
1727 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: KADMM*S' ", e_dummy
1728 : END IF
1729 : END IF
1730 :
1731 162 : ALLOCATE (matrix_ks_aux(nspins))
1732 90 : DO ispin = 1, nspins
1733 54 : NULLIFY (matrix_ks_aux(ispin)%matrix)
1734 54 : ALLOCATE (matrix_ks_aux(ispin)%matrix)
1735 54 : CALL dbcsr_copy(matrix_ks_aux(ispin)%matrix, p_env%kpp1_admm(ispin)%matrix)
1736 90 : CALL dbcsr_set(matrix_ks_aux(ispin)%matrix, 0.0_dp)
1737 : END DO
1738 :
1739 : ! Calculate kernel
1740 36 : CALL tddft_hfx_matrix(matrix_ks_aux, rho_ao_aux, qs_env, .FALSE.)
1741 :
1742 36 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
1743 24 : CALL get_qs_env(qs_env, ks_env=ks_env)
1744 24 : CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list_aux_fit)
1745 :
1746 60 : DO ispin = 1, nspins
1747 60 : CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
1748 : END DO
1749 :
1750 24 : IF (use_virial) THEN
1751 8 : CALL qs_rho_get(p_env%rho1_admm, rho_r=rho_mp2_r_aux)
1752 8 : e_xc = 0.0_dp
1753 20 : DO ispin = 1, nspins
1754 20 : e_xc = e_xc + pw_integral_ab(rho_mp2_r_aux(ispin), vadmm_rspace(ispin))
1755 : END DO
1756 :
1757 8 : e_xc = -e_xc/vadmm_rspace(1)%pw_grid%dvol/REAL(para_env%num_pe, dp)
1758 :
1759 : ! Update the virial
1760 32 : DO alpha = 1, 3
1761 24 : virial%pv_exc(alpha, alpha) = virial%pv_exc(alpha, alpha) + e_xc
1762 32 : virial%pv_virial(alpha, alpha) = virial%pv_virial(alpha, alpha) + e_xc
1763 : END DO
1764 : IF (debug_forces) THEN
1765 8 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: P1*VADMM ", e_xc
1766 : END IF
1767 : END IF
1768 :
1769 120 : IF (use_virial) h_stress = virial%pv_virial
1770 : IF (debug_forces) THEN
1771 96 : deb = force(1)%rho_elec(1:3, 1)
1772 24 : IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1773 : END IF
1774 60 : DO ispin = 1, nspins
1775 36 : IF (do_exx) THEN
1776 : CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=matrix_ks_aux(ispin), qs_env=qs_env, &
1777 : calculate_forces=.TRUE., basis_type="AUX_FIT", &
1778 0 : task_list_external=task_list_aux_fit, pmat=matrix_p_mp2_admm(ispin))
1779 : ELSE
1780 : CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=matrix_ks_aux(ispin), qs_env=qs_env, &
1781 : calculate_forces=.TRUE., basis_type="AUX_FIT", &
1782 36 : task_list_external=task_list_aux_fit, pmat=rho_ao_aux(ispin))
1783 : END IF
1784 60 : CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
1785 : END DO
1786 120 : IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - h_stress)
1787 24 : DEALLOCATE (vadmm_rspace)
1788 : IF (debug_forces) THEN
1789 96 : deb(:) = force(1)%rho_elec(:, 1) - deb
1790 24 : CALL para_env%sum(deb)
1791 24 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*VADMM' ", deb
1792 24 : IF (use_virial) THEN
1793 8 : e_dummy = third_tr(virial%pv_virial) - e_dummy
1794 8 : CALL para_env%sum(e_dummy)
1795 8 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: VADMM' ", e_dummy
1796 : END IF
1797 : END IF
1798 :
1799 60 : DO ispin = 1, nspins
1800 60 : CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
1801 : END DO
1802 :
1803 : END IF
1804 :
1805 90 : DO ispin = 1, nspins
1806 90 : CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
1807 : END DO
1808 :
1809 : IF (debug_forces) THEN
1810 144 : deb = force(1)%overlap_admm(1:3, 1)
1811 36 : IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1812 : END IF
1813 : ! Add the second half of the projector deriatives contracting the first order density matrix with the fockian in the auxiliary basis
1814 36 : IF (do_exx) THEN
1815 0 : CALL admm_projection_derivative(qs_env, matrix_ks_aux, matrix_p_mp2)
1816 : ELSE
1817 36 : CALL admm_projection_derivative(qs_env, matrix_ks_aux, rho_ao)
1818 : END IF
1819 : IF (debug_forces) THEN
1820 144 : deb(:) = force(1)%overlap_admm(:, 1) - deb
1821 36 : CALL para_env%sum(deb)
1822 36 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*VADMM*dS'", deb
1823 36 : IF (use_virial) THEN
1824 12 : e_dummy = third_tr(virial%pv_virial) - e_dummy
1825 12 : CALL para_env%sum(e_dummy)
1826 12 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: VADMM*S' ", e_dummy
1827 : END IF
1828 : END IF
1829 :
1830 90 : DO ispin = 1, nspins
1831 90 : CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, -1.0_dp)
1832 : END DO
1833 :
1834 90 : DO ispin = 1, nspins
1835 54 : CALL dbcsr_release(matrix_ks_aux(ispin)%matrix)
1836 90 : DEALLOCATE (matrix_ks_aux(ispin)%matrix)
1837 : END DO
1838 36 : DEALLOCATE (matrix_ks_aux)
1839 : END IF
1840 : END IF
1841 :
1842 272 : CALL dbcsr_scale(p_env%w1(1)%matrix, -1.0_dp)
1843 :
1844 : ! Finish matrix_w_mp2 with occ-occ block
1845 630 : DO ispin = 1, nspins
1846 358 : CALL get_mo_set(mo_set=mos(ispin), homo=nocc, nmo=alpha)
1847 : CALL calculate_whz_matrix(mos(ispin)%mo_coeff, p_env%kpp1(ispin)%matrix, &
1848 630 : p_env%w1(1)%matrix, 1.0_dp, nocc)
1849 : END DO
1850 :
1851 272 : IF (debug_forces .AND. use_virial) e_dummy = third_tr(virial%pv_virial)
1852 :
1853 272 : NULLIFY (scrm)
1854 : CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
1855 : matrix_name="OVERLAP MATRIX", &
1856 : basis_type_a="ORB", basis_type_b="ORB", &
1857 : sab_nl=sab_orb, calculate_forces=.TRUE., &
1858 272 : matrix_p=p_env%w1(1)%matrix)
1859 272 : CALL dbcsr_deallocate_matrix_set(scrm)
1860 :
1861 : IF (debug_forces) THEN
1862 1088 : deb = force(1)%overlap(1:3, 1)
1863 272 : CALL para_env%sum(deb)
1864 272 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: W*dS ", deb
1865 272 : IF (use_virial) THEN
1866 50 : e_dummy = third_tr(virial%pv_virial) - e_dummy
1867 50 : CALL para_env%sum(e_dummy)
1868 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: S ", e_dummy
1869 : END IF
1870 : END IF
1871 :
1872 272 : CALL auxbas_pw_pool%give_back_pw(pot_r)
1873 272 : CALL auxbas_pw_pool%give_back_pw(pot_g)
1874 272 : CALL auxbas_pw_pool%give_back_pw(rho_tot_g)
1875 :
1876 : ! Release linres stuff
1877 272 : CALL p_env_release(p_env)
1878 272 : DEALLOCATE (p_env)
1879 272 : NULLIFY (qs_env%mp2_env%ri_grad%p_env)
1880 :
1881 272 : CALL sum_qs_force(force, qs_env%mp2_env%ri_grad%mp2_force)
1882 272 : CALL deallocate_qs_force(qs_env%mp2_env%ri_grad%mp2_force)
1883 :
1884 272 : IF (use_virial) THEN
1885 1250 : virial%pv_mp2 = qs_env%mp2_env%ri_grad%mp2_virial
1886 1250 : virial%pv_virial = virial%pv_virial + qs_env%mp2_env%ri_grad%mp2_virial
1887 : IF (debug_forces) THEN
1888 50 : e_dummy = third_tr(virial%pv_mp2)
1889 50 : CALL para_env%sum(e_dummy)
1890 50 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: MP2nonsep ", e_dummy
1891 : END IF
1892 : END IF
1893 : ! Rewind the change from the beginning
1894 272 : IF (use_virial) virial%pv_calculate = .FALSE.
1895 :
1896 : ! Add the dispersion forces and virials
1897 272 : CALL get_qs_env(qs_env, dispersion_env=dispersion_env)
1898 272 : CALL calculate_dispersion_pairpot(qs_env, dispersion_env, e_dummy, .TRUE.)
1899 :
1900 : CALL cp_print_key_finished_output(iounit, logger, input, &
1901 272 : "DFT%XC%WF_CORRELATION%PRINT")
1902 :
1903 272 : CALL timestop(handle)
1904 :
1905 544 : END SUBROUTINE update_mp2_forces
1906 :
1907 : ! **************************************************************************************************
1908 : !> \brief Calculates the third of the trace of a 3x3 matrix, for debugging purposes
1909 : !> \param matrix ...
1910 : !> \return ...
1911 : ! **************************************************************************************************
1912 665 : PURE FUNCTION third_tr(matrix)
1913 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN) :: matrix
1914 : REAL(KIND=dp) :: third_tr
1915 :
1916 665 : third_tr = (matrix(1, 1) + matrix(2, 2) + matrix(3, 3))/3.0_dp
1917 :
1918 665 : END FUNCTION third_tr
1919 :
1920 : END MODULE mp2_cphf
|
