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