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