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