(git:e7e05ae)
mp2_cphf.F
Go to the documentation of this file.
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  SUBROUTINE solve_z_vector_eq(qs_env, mp2_env, para_env, dft_control, &
160  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  INTEGER, ALLOCATABLE, DIMENSION(:) :: virtual
176  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  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: l_jb, mo_coeff_o, mo_coeff_v, p_ia, &
184  p_mo, w_mo
185  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p_mp2, &
186  matrix_p_mp2_admm, matrix_s, p_mu_nu, &
187  rho1_ao, rho_ao, rho_ao_aux_fit
188  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  CALL timeset(routinen, handle)
198 
199  ! start collecting stuff
200  dimen = nmo
201  NULLIFY (input, matrix_s, blacs_env, rho, &
202  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  rho=rho)
212 
213  CALL qs_rho_get(rho, rho_ao=rho_ao)
214 
215  ! Get number of relevant spin states
216  nspins = dft_control%nspins
217  alpha_beta = (nspins == 2)
218 
219  CALL move_alloc(mp2_env%ri_grad%P_mo, p_mo)
220  CALL move_alloc(mp2_env%ri_grad%W_mo, w_mo)
221  CALL move_alloc(mp2_env%ri_grad%L_jb, l_jb)
222 
223  ALLOCATE (virtual(nspins))
224  virtual(:) = dimen - homo(:)
225 
226  NULLIFY (p_mu_nu)
227  CALL dbcsr_allocate_matrix_set(p_mu_nu, nspins)
228  DO ispin = 1, nspins
229  ALLOCATE (p_mu_nu(ispin)%matrix)
230  CALL dbcsr_copy(p_mu_nu(ispin)%matrix, rho_ao(1)%matrix, name="P_mu_nu")
231  CALL dbcsr_set(p_mu_nu(ispin)%matrix, 0.0_dp)
232  END DO
233 
234  NULLIFY (fm_struct_tmp)
235  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
236  nrow_global=dimen, ncol_global=dimen)
237  CALL cp_fm_create(fm_g_mu_nu, fm_struct_tmp, name="G_mu_nu")
238  CALL cp_fm_create(fm_back, fm_struct_tmp, name="fm_back")
239  CALL cp_fm_struct_release(fm_struct_tmp)
240  CALL cp_fm_set_all(fm_g_mu_nu, 0.0_dp)
241  CALL cp_fm_set_all(fm_back, 0.0_dp)
242 
243  ALLOCATE (mo_coeff_o(nspins), mo_coeff_v(nspins))
244  DO ispin = 1, nspins
245  NULLIFY (fm_struct_tmp)
246  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
247  nrow_global=dimen, ncol_global=homo(ispin))
248  CALL cp_fm_create(mo_coeff_o(ispin), fm_struct_tmp, name="mo_coeff_o")
249  CALL cp_fm_struct_release(fm_struct_tmp)
250  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  t_firstrow=1, t_firstcol=1)
255 
256  NULLIFY (fm_struct_tmp)
257  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
258  nrow_global=dimen, ncol_global=virtual(ispin))
259  CALL cp_fm_create(mo_coeff_v(ispin), fm_struct_tmp, name="mo_coeff_v")
260  CALL cp_fm_struct_release(fm_struct_tmp)
261  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  t_firstrow=1, t_firstcol=1)
266  END DO
267 
268  ! hfx section
269  NULLIFY (hfx_sections)
270  hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
271  CALL section_vals_get(hfx_sections, explicit=do_hfx, n_repetition=n_rep_hf)
272  IF (do_hfx) THEN
273  ! here we check if we have to reallocate the HFX container
274  IF (mp2_env%ri_grad%free_hfx_buffer .AND. (.NOT. qs_env%x_data(1, 1)%do_hfx_ri)) THEN
275  CALL timeset(routinen//"_alloc_hfx", handle2)
276  n_threads = 1
277 !$ n_threads = omp_get_max_threads()
278 
279  DO irep = 1, n_rep_hf
280  DO i_thread = 0, n_threads - 1
281  actual_x_data => qs_env%x_data(irep, i_thread + 1)
282 
283  do_dynamic_load_balancing = .true.
284  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  my_bin_size = SIZE(actual_x_data%distribution_energy)
288  ELSE
289  my_bin_size = 1
290  END IF
291 
292  IF (.NOT. actual_x_data%memory_parameter%do_all_on_the_fly) THEN
293  CALL alloc_containers(actual_x_data%store_ints, my_bin_size)
294 
295  DO bin = 1, my_bin_size
296  maxval_container => actual_x_data%store_ints%maxval_container(bin)
297  integral_containers => actual_x_data%store_ints%integral_containers(:, bin)
298  CALL hfx_init_container(maxval_container, actual_x_data%memory_parameter%actual_memory_usage, .false.)
299  DO i = 1, 64
300  CALL hfx_init_container(integral_containers(i), &
301  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  CALL timestop(handle2)
308  END IF
309 
310  ! set up parameters for P_screening
311  restore_p_screen = qs_env%x_data(1, 1)%screening_parameter%do_initial_p_screening
312  IF (qs_env%x_data(1, 1)%screening_parameter%do_initial_p_screening) THEN
313  IF (mp2_env%ri_grad%free_hfx_buffer) THEN
314  mp2_env%p_screen = .false.
315  ELSE
316  mp2_env%p_screen = .true.
317  END IF
318  END IF
319  END IF
320 
321  ! Add exx part for RPA
322  do_exx = .false.
323  IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
324  hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
325  CALL section_vals_get(hfx_section, explicit=do_exx)
326  END IF
327  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  reuse_hfx=mp2_env%ri_rpa%reuse_hfx)
336 
337  focc = 1.0_dp
338  IF (nspins == 1) focc = 2.0_dp
339  !focc = 0.0_dp
340  DO ispin = 1, nspins
341  CALL dbcsr_add(p_mu_nu(ispin)%matrix, matrix_ks(ispin)%matrix, 1.0_dp, -1.0_dp)
342  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  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  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  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  ALLOCATE (linres_control)
355  linres_control%do_kernel = .true.
356  linres_control%lr_triplet = .false.
357  linres_control%linres_restart = .false.
358  linres_control%max_iter = mp2_env%ri_grad%cphf_max_num_iter
359  linres_control%eps = mp2_env%ri_grad%cphf_eps_conv
360  linres_control%eps_filter = mp2_env%mp2_gpw%eps_filter
361  linres_control%restart_every = 50
362  linres_control%preconditioner_type = ot_precond_full_all
363  linres_control%energy_gap = 0.02_dp
364 
365  NULLIFY (p_env)
366  ALLOCATE (p_env)
367  CALL p_env_create(p_env, qs_env, p1_option=p_mu_nu, orthogonal_orbitals=.true., linres_control=linres_control)
368  CALL set_qs_env(qs_env, linres_control=linres_control)
369  CALL p_env_psi0_changed(p_env, qs_env)
370  p_env%new_preconditioner = .true.
371  CALL p_env_check_i_alloc(p_env, qs_env)
372  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  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  .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  ALLOCATE (p_ia(nspins))
393  DO ispin = 1, nspins
394  NULLIFY (fm_struct_tmp)
395  CALL cp_fm_struct_create(fm_struct_tmp, para_env=para_env, context=blacs_env, &
396  nrow_global=homo(ispin), ncol_global=virtual(ispin))
397  CALL cp_fm_create(p_ia(ispin), fm_struct_tmp, name="P_ia")
398  CALL cp_fm_struct_release(fm_struct_tmp)
399  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  l_jb, fm_g_mu_nu, fm_back, p_ia)
405 
406  ! release fm stuff
407  CALL cp_fm_release(fm_g_mu_nu)
408  CALL cp_fm_release(l_jb)
409  CALL cp_fm_release(mo_coeff_o)
410  CALL cp_fm_release(mo_coeff_v)
411 
412  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  t_firstrow=1, t_firstcol=homo(ispin) + 1)
418  ! transpose P_MO matrix (easy way to symmetrize)
419  CALL cp_fm_set_all(fm_back, 0.0_dp)
420  ! P_mo now is ready
421  CALL cp_fm_upper_to_full(matrix=p_mo(ispin), work=fm_back)
422  END DO
423  CALL cp_fm_release(p_ia)
424 
425  ! do the final update to the energy weighted matrix W_MO
426  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  col_indices=col_indices)
432  DO jjb = 1, ncol_local
433  j_global = col_indices(jjb)
434  IF (j_global <= homo(ispin)) THEN
435  DO iib = 1, nrow_local
436  i_global = row_indices(iib)
437  w_mo(ispin)%local_data(iib, jjb) = w_mo(ispin)%local_data(iib, jjb) &
438  - p_mo(ispin)%local_data(iib, jjb)*eigenval(j_global, ispin)
439  IF (i_global == j_global .AND. nspins == 1) w_mo(ispin)%local_data(iib, jjb) = &
440  w_mo(ispin)%local_data(iib, jjb) - 2.0_dp*eigenval(j_global, ispin)
441  IF (i_global == j_global .AND. nspins == 2) w_mo(ispin)%local_data(iib, jjb) = &
442  w_mo(ispin)%local_data(iib, jjb) - eigenval(j_global, ispin)
443  END DO
444  ELSE
445  DO iib = 1, nrow_local
446  i_global = row_indices(iib)
447  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  - 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  - 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  NULLIFY (p_env%w1)
463  CALL dbcsr_allocate_matrix_set(p_env%w1, 1)
464  ALLOCATE (p_env%w1(1)%matrix)
465  CALL dbcsr_copy(p_env%w1(1)%matrix, matrix_s(1)%matrix, &
466  name="W MATRIX MP2")
467  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  DO ispin = 1, nspins
471  CALL parallel_gemm('N', 'N', dimen, dimen, dimen, 1.0_dp, &
472  mo_coeff(ispin), w_mo(ispin), 0.0_dp, fm_back)
473  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  CALL cp_fm_release(w_mo)
476 
477  CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
478 
479  DO ispin = 1, nspins
480  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  mo_coeff(ispin), p_mo(ispin), 0.0_dp, fm_back)
484  CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff(ispin), dimen, 1.0_dp, .true.)
485 
486  CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix)
487  END DO
488  CALL cp_fm_release(p_mo)
489  CALL cp_fm_release(fm_back)
490 
491  CALL p_env_update_rho(p_env, qs_env)
492 
493  ! create mp2 DBCSR density
494  CALL dbcsr_allocate_matrix_set(matrix_p_mp2, nspins)
495  DO ispin = 1, nspins
496  ALLOCATE (matrix_p_mp2(ispin)%matrix)
497  CALL dbcsr_copy(matrix_p_mp2(ispin)%matrix, p_env%p1(ispin)%matrix, &
498  name="P MATRIX MP2")
499  END DO
500 
501  IF (dft_control%do_admm) THEN
502  CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux_fit)
503  CALL qs_rho_get(rho_aux_fit, rho_ao=rho_ao_aux_fit)
504 
505  ! create mp2 DBCSR density in auxiliary basis
506  CALL dbcsr_allocate_matrix_set(matrix_p_mp2_admm, nspins)
507  DO ispin = 1, nspins
508  ALLOCATE (matrix_p_mp2_admm(ispin)%matrix)
509  CALL dbcsr_copy(matrix_p_mp2_admm(ispin)%matrix, p_env%p1_admm(ispin)%matrix, &
510  name="P MATRIX MP2 ADMM")
511  END DO
512  END IF
513 
514  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  mp2_env%not_last_hfx = .false.
518  mp2_env%p_screen = restore_p_screen
519 
520  CALL timestop(handle)
521 
522  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  SUBROUTINE cphf_like_update(qs_env, homo, virtual, dimen, nspins, &
549  mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
550  fm_mo, fm_ao, fm_back, transf_type_in, &
551  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  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
570  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho1_ao
571 
572  CALL timeset(routinen, handle)
573 
574  CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
575 
576  ! Determine the first-order density matrices in AO basis
577  DO ispin = 1, nspins
578  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  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  mo_coeff_o(ispin), mat_in, 0.0_dp, fm_back)
588  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  b_first_row=homo(ispin) + 1)
594  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  mo_coeff_o(ispin), mat_in, 0.0_dp, fm_back)
600  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  mo_coeff_v(ispin), mat_in, 0.0_dp, fm_back)
604  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  CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix)
612  END DO
613 
614  CALL p_env_update_rho(p_env, qs_env)
615 
616  IF (transf_type_in == 3) THEN
617  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  col_indices=col_indices)
625  DO jjb = 1, ncol_local
626  j_global = col_indices(jjb)
627  DO iib = 1, nrow_local
628  i_global = row_indices(iib)
629  fm_mo_out(ispin)%local_data(iib, jjb) = fm_mo(ispin)%local_data(iib, jjb)* &
630  (eigenval(j_global + homo(ispin), ispin) - eigenval(i_global, ispin))
631  END DO
632  END DO
633  END DO
634  END IF
635 
636  CALL apply_2nd_order_kernel(qs_env, p_env, recalc_hfx_integrals)
637 
638  DO ispin = 1, nspins
639  ! copy back to fm
640  CALL cp_fm_set_all(fm_ao, 0.0_dp)
641  CALL copy_dbcsr_to_fm(matrix=p_env%kpp1(ispin)%matrix, fm=fm_ao)
642  CALL cp_fm_set_all(fm_back, 0.0_dp)
643  CALL cp_fm_upper_to_full(fm_ao, fm_back)
644 
645  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  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  fm_back, mo_coeff_v(ispin), 1.0_dp, mat_out)
654  END associate
655  END DO
656 
657  CALL timestop(handle)
658 
659  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  SUBROUTINE solve_z_vector_eq_low(qs_env, mp2_env, homo, virtual, dimen, unit_nr, nspins, &
682  mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
683  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  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
701  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: precond
702 
703  CALL timeset(routinen, handle)
704 
705  ! Pople method
706  ! change sign to L_jb
707  DO ispin = 1, nspins
708  l_jb(ispin)%local_data(:, :) = -l_jb(ispin)%local_data(:, :)
709  END DO
710 
711  ! create fm structure
712  ALLOCATE (precond(nspins))
713  DO ispin = 1, nspins
714  ! create preconditioner (for now only orbital energy differences)
715  CALL cp_fm_create(precond(ispin), p_ia(ispin)%matrix_struct, name="precond")
716  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  col_indices=col_indices)
722  DO jjb = 1, ncol_local
723  j_global = col_indices(jjb)
724  DO iib = 1, nrow_local
725  i_global = row_indices(iib)
726  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  DO ispin = 1, nspins
732  precond(ispin)%local_data(:, :) = 1.0_dp/precond(ispin)%local_data(:, :)
733  END DO
734 
735  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  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  l_jb, fm_g_mu_nu, fm_back, p_ia, precond)
744  CASE DEFAULT
745  cpabort("Unknown solver")
746  END SELECT
747 
748  CALL cp_fm_release(precond)
749 
750  CALL timestop(handle)
751 
752  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  SUBROUTINE solve_z_vector_pople(qs_env, mp2_env, homo, virtual, dimen, unit_nr, nspins, &
774  mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
775  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: proj_bi_xj, temp_vals, x_norm, xi_b
795  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a_small, b_small, xi_axi
796  TYPE(cp_fm_struct_p_type), ALLOCATABLE, &
797  DIMENSION(:) :: fm_struct_tmp
798  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: b_i, residual
799  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: ax, xn
800 
801  CALL timeset(routinen, handle)
802 
803  eps_conv = mp2_env%ri_grad%cphf_eps_conv
804 
805  IF (sqrt(accurate_dot_product_spin(l_jb, l_jb)) >= eps_conv) THEN
806 
807  max_num_iter = mp2_env%ri_grad%cphf_max_num_iter
808  scale_cphf = mp2_env%ri_grad%scale_step_size
809 
810  ! set the transformation type (equal for all methods all updates)
811  transf_type_in = 3
812 
813  ! set convergence flag
814  converged = .false.
815 
816  ALLOCATE (fm_struct_tmp(nspins), b_i(nspins), residual(nspins))
817  DO ispin = 1, nspins
818  fm_struct_tmp(ispin)%struct => p_ia(ispin)%matrix_struct
819 
820  CALL cp_fm_create(b_i(ispin), fm_struct_tmp(ispin)%struct, name="b_i")
821  CALL cp_fm_set_all(b_i(ispin), 0.0_dp)
822  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  CALL cp_fm_create(residual(ispin), fm_struct_tmp(ispin)%struct, name="residual")
827  CALL cp_fm_set_all(residual(ispin), 0.0_dp)
828  END DO
829 
830  IF (unit_nr > 0) THEN
831  WRITE (unit_nr, *)
832  WRITE (unit_nr, '(T3,A)') "MP2_CPHF| Iterative solution of Z-Vector equations (Pople's method)"
833  WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', eps_conv
834  WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', max_num_iter
835  WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Scaling of initial guess: ', scale_cphf
836  WRITE (unit_nr, '(T4,A)') repeat("-", 40)
837  WRITE (unit_nr, '(T4,A,T15,A,T33,A)') 'Step', 'Time', 'Convergence'
838  WRITE (unit_nr, '(T4,A)') repeat("-", 40)
839  END IF
840 
841  ALLOCATE (xn(nspins, max_num_iter))
842  ALLOCATE (ax(nspins, max_num_iter))
843  ALLOCATE (x_norm(max_num_iter))
844  ALLOCATE (xi_b(max_num_iter))
845  ALLOCATE (xi_axi(max_num_iter, 0:max_num_iter))
846  x_norm = 0.0_dp
847  xi_b = 0.0_dp
848  xi_axi = 0.0_dp
849 
850  cycle_counter = 0
851  DO iiter = 1, max_num_iter
852  cycle_counter = cycle_counter + 1
853 
854  t1 = m_walltime()
855 
856  ! create and update x_i (orthogonalization with previous vectors)
857  DO ispin = 1, nspins
858  CALL cp_fm_create(xn(ispin, iiter), fm_struct_tmp(ispin)%struct, name="xi")
859  CALL cp_fm_set_all(xn(ispin, iiter), 0.0_dp)
860  END DO
861 
862  ALLOCATE (proj_bi_xj(iiter - 1))
863  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  DO iib = 1, iiter - 1
867  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  DO ispin = 1, nspins
872  xn(ispin, iiter)%local_data(:, :) = scale_cphf*b_i(ispin)%local_data(:, :)
873  DO iib = 1, iiter - 1
874  xn(ispin, iiter)%local_data(:, :) = xn(ispin, iiter)%local_data(:, :) - &
875  xn(ispin, iib)%local_data(:, :)*proj_bi_xj(iib)
876  END DO
877  END DO
878  DEALLOCATE (proj_bi_xj)
879 
880  ! create Ax(iiter) that will store the matrix vector product for this cycle
881  DO ispin = 1, nspins
882  CALL cp_fm_create(ax(ispin, iiter), fm_struct_tmp(ispin)%struct, name="Ai")
883  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  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  ALLOCATE (temp_vals(iiter + 2))
900  temp_vals = 0.0_dp
901  ! <Ax_i|x_j>
902  DO iib = 1, iiter
903  temp_vals(iib) = temp_vals(iib) + accurate_dot_product_spin(ax(:, iiter), xn(:, iib))
904  END DO
905  ! <x_i|b>
906  temp_vals(iiter + 1) = temp_vals(iiter + 1) + accurate_dot_product_spin(xn(:, iiter), l_jb)
907  ! norm
908  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  xi_axi(iiter, 1:iiter) = temp_vals(1:iiter)
911  xi_axi(1:iiter, iiter) = temp_vals(1:iiter)
912  xi_b(iiter) = temp_vals(iiter + 1)
913  x_norm(iiter) = temp_vals(iiter + 2)
914  DEALLOCATE (temp_vals)
915 
916  ! solve reduced system
917  IF (ALLOCATED(a_small)) DEALLOCATE (a_small)
918  IF (ALLOCATED(b_small)) DEALLOCATE (b_small)
919  ALLOCATE (a_small(iiter, iiter))
920  ALLOCATE (b_small(iiter, 1))
921  a_small(1:iiter, 1:iiter) = xi_axi(1:iiter, 1:iiter)
922  b_small(1:iiter, 1) = xi_b(1:iiter)
923 
924  CALL solve_system(matrix=a_small, mysize=iiter, eigenvectors=b_small)
925 
926  ! check for convergence
927  DO ispin = 1, nspins
928  CALL cp_fm_set_all(residual(ispin), 0.0_dp)
929  DO iib = 1, iiter
930  residual(ispin)%local_data(:, :) = &
931  residual(ispin)%local_data(:, :) + &
932  b_small(iib, 1)*ax(ispin, iib)%local_data(:, :)
933  END DO
934 
935  residual(ispin)%local_data(:, :) = &
936  residual(ispin)%local_data(:, :) - &
937  l_jb(ispin)%local_data(:, :)
938  END DO
939 
940  conv = sqrt(accurate_dot_product_spin(residual, residual))
941 
942  t2 = m_walltime()
943 
944  IF (unit_nr > 0) THEN
945  WRITE (unit_nr, '(T3,I5,T13,F6.1,11X,F14.8)') iiter, t2 - t1, conv
946  CALL m_flush(unit_nr)
947  END IF
948 
949  IF (conv <= eps_conv) THEN
950  converged = .true.
951  EXIT
952  END IF
953 
954  ! update b_i for the next round
955  DO ispin = 1, nspins
956  b_i(ispin)%local_data(:, :) = b_i(ispin)%local_data(:, :) &
957  + precond(ispin)%local_data(:, :) &
958  *ax(ispin, iiter)%local_data(:, :)
959  END DO
960 
961  scale_cphf = 1.0_dp
962 
963  END DO
964 
965  IF (unit_nr > 0) THEN
966  WRITE (unit_nr, '(T4,A)') repeat("-", 40)
967  IF (converged) THEN
968  WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations converged in', cycle_counter, ' steps'
969  ELSE
970  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  DO iiter = 1, cycle_counter
976  DO ispin = 1, nspins
977  p_ia(ispin)%local_data(:, :) = p_ia(ispin)%local_data(:, :) + &
978  b_small(iiter, 1)*xn(ispin, iiter)%local_data(:, :)
979  END DO
980  END DO
981 
982  ! Release arrays
983  DEALLOCATE (x_norm)
984  DEALLOCATE (xi_b)
985  DEALLOCATE (xi_axi)
986 
987  CALL cp_fm_release(b_i)
988  CALL cp_fm_release(residual)
989  CALL cp_fm_release(ax)
990  CALL cp_fm_release(xn)
991 
992  ELSE
993  IF (unit_nr > 0) THEN
994  WRITE (unit_nr, '(T4,A)') repeat("-", 40)
995  WRITE (unit_nr, '(T3,A)') 'Residual smaller than EPS_CONV. Skip solution of Z-vector equation.'
996  END IF
997  END IF
998 
999  CALL timestop(handle)
1000 
1001  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  SUBROUTINE solve_z_vector_cg(qs_env, mp2_env, homo, virtual, dimen, unit_nr, nspins, &
1023  mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
1024  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  DIMENSION(:) :: fm_struct_tmp
1048  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: a_dot_search_vector, diff_search_vector, &
1049  residual, search_vector
1050 
1051  CALL timeset(routinen, handle)
1052 
1053  max_num_iter = mp2_env%ri_grad%cphf_max_num_iter
1054  eps_conv = mp2_env%ri_grad%cphf_eps_conv
1055  z_solver_method = mp2_env%ri_grad%z_solver_method
1056  restart_every = mp2_env%ri_grad%cphf_restart
1057  scale_step_size = mp2_env%ri_grad%scale_step_size
1058  transf_type_in = 3
1059 
1060  IF (unit_nr > 0) THEN
1061  WRITE (unit_nr, *)
1062  SELECT CASE (z_solver_method)
1063  CASE (z_solver_cg)
1064  IF (mp2_env%ri_grad%polak_ribiere) THEN
1065  WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (CG with Polak-Ribiere step)'
1066  ELSE
1067  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  WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (Richardson method)'
1071  CASE (z_solver_sd)
1072  WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (Steepest Descent method)'
1073  CASE DEFAULT
1074  cpabort("Unknown solver")
1075  END SELECT
1076  WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', eps_conv
1077  WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', max_num_iter
1078  WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Number of steps for restart: ', restart_every
1079  WRITE (unit_nr, '(T3, A)') 'MP2_CPHF| Restart after no decrease'
1080  WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Scaling factor of each step: ', scale_step_size
1081  WRITE (unit_nr, '(T4,A)') repeat("-", 40)
1082  WRITE (unit_nr, '(T4,A,T13,A,T28,A,T43,A)') 'Step', 'Restart', 'Time', 'Convergence'
1083  WRITE (unit_nr, '(T4,A)') repeat("-", 40)
1084  END IF
1085 
1086  ALLOCATE (fm_struct_tmp(nspins), residual(nspins), diff_search_vector(nspins), &
1087  search_vector(nspins), a_dot_search_vector(nspins))
1088  DO ispin = 1, nspins
1089  fm_struct_tmp(ispin)%struct => p_ia(ispin)%matrix_struct
1090 
1091  CALL cp_fm_create(residual(ispin), fm_struct_tmp(ispin)%struct, name="residual")
1092  CALL cp_fm_set_all(residual(ispin), 0.0_dp)
1093 
1094  CALL cp_fm_create(diff_search_vector(ispin), fm_struct_tmp(ispin)%struct, name="difference search vector")
1095  CALL cp_fm_set_all(diff_search_vector(ispin), 0.0_dp)
1096 
1097  CALL cp_fm_create(search_vector(ispin), fm_struct_tmp(ispin)%struct, name="search vector")
1098  CALL cp_fm_set_all(search_vector(ispin), 0.0_dp)
1099 
1100  CALL cp_fm_create(a_dot_search_vector(ispin), fm_struct_tmp(ispin)%struct, name="A times search vector")
1101  CALL cp_fm_set_all(a_dot_search_vector(ispin), 0.0_dp)
1102  END DO
1103 
1104  converged = .false.
1105  cycles_passed = max_num_iter
1106  ! By that, we enforce the setup of the matrices
1107  do_restart = .true.
1108 
1109  t1 = m_walltime()
1110 
1111  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  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  IF (.NOT. mp2_env%ri_grad%recalc_residual .OR. (iiter == 1)) THEN
1118  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  residual)
1124  ELSE
1125  do_restart = .false.
1126 
1127  DO ispin = 1, nspins
1128  CALL cp_fm_set_all(residual(ispin), 0.0_dp)
1129  END DO
1130  END IF
1131 
1132  DO ispin = 1, nspins
1133  residual(ispin)%local_data(:, :) = l_jb(ispin)%local_data(:, :) &
1134  - residual(ispin)%local_data(:, :)
1135  END DO
1136  END IF
1137 
1138  DO ispin = 1, nspins
1139  diff_search_vector(ispin)%local_data(:, :) = &
1140  precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1141  search_vector(ispin)%local_data(:, :) = diff_search_vector(ispin)%local_data(:, :)
1142  END DO
1143 
1144  restart_counter = 1
1145  END IF
1146 
1147  norm_residual_old = sqrt(accurate_dot_product_spin(residual, residual))
1148 
1149  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  a_dot_search_vector)
1156 
1157  IF (z_solver_method /= z_solver_richardson) THEN
1158  search_vec_dot_a_search_vec = accurate_dot_product_spin(search_vector, a_dot_search_vector)
1159 
1160  IF (z_solver_method == z_solver_cg) THEN
1161  scale_result = residual_dot_diff_search_vec_old/search_vec_dot_a_search_vec
1162  ELSE
1163  residual_dot_search_vec = accurate_dot_product_spin(residual, search_vector)
1164  scale_result = residual_dot_search_vec/search_vec_dot_a_search_vec
1165  END IF
1166 
1167  scale_result = scale_result*scale_step_size
1168 
1169  ELSE
1170 
1171  scale_result = scale_step_size
1172 
1173  END IF
1174 
1175  DO ispin = 1, nspins
1176  p_ia(ispin)%local_data(:, :) = p_ia(ispin)%local_data(:, :) &
1177  + scale_result*search_vector(ispin)%local_data(:, :)
1178  END DO
1179 
1180  IF (.NOT. mp2_env%ri_grad%recalc_residual) THEN
1181 
1182  DO ispin = 1, nspins
1183  residual(ispin)%local_data(:, :) = residual(ispin)%local_data(:, :) &
1184  - 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  residual)
1192 
1193  DO ispin = 1, nspins
1194  residual(ispin)%local_data(:, :) = l_jb(ispin)%local_data(:, :) - residual(ispin)%local_data(:, :)
1195  END DO
1196  END IF
1197 
1198  norm_residual = sqrt(accurate_dot_product_spin(residual, residual))
1199 
1200  t2 = m_walltime()
1201 
1202  IF (unit_nr > 0) THEN
1203  WRITE (unit_nr, '(T3,I4,T16,L1,T26,F6.1,8X,F14.8)') iiter, do_restart, t2 - t1, norm_residual
1204  CALL m_flush(unit_nr)
1205  END IF
1206 
1207  IF (norm_residual <= eps_conv) THEN
1208  converged = .true.
1209  cycles_passed = iiter
1210  EXIT
1211  END IF
1212 
1213  t1 = m_walltime()
1214 
1215  IF (z_solver_method == z_solver_richardson) THEN
1216  DO ispin = 1, nspins
1217  search_vector(ispin)%local_data(:, :) = &
1218  scale_step_size*precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1219  END DO
1220  ELSE IF (z_solver_method == z_solver_sd) THEN
1221  DO ispin = 1, nspins
1222  search_vector(ispin)%local_data(:, :) = &
1223  precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1224  END DO
1225  ELSE
1226  IF (mp2_env%ri_grad%polak_ribiere) &
1227  residual_new_dot_diff_search_vec_old = accurate_dot_product_spin(residual, diff_search_vector)
1228 
1229  DO ispin = 1, nspins
1230  diff_search_vector(ispin)%local_data(:, :) = &
1231  precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1232  END DO
1233 
1234  residual_dot_diff_search_vec_new = accurate_dot_product_spin(residual, diff_search_vector)
1235 
1236  scale_search = residual_dot_diff_search_vec_new/residual_dot_diff_search_vec_old
1237  IF (mp2_env%ri_grad%polak_ribiere) scale_search = &
1238  scale_search - residual_new_dot_diff_search_vec_old/residual_dot_diff_search_vec_old
1239 
1240  DO ispin = 1, nspins
1241  search_vector(ispin)%local_data(:, :) = scale_search*search_vector(ispin)%local_data(:, :) &
1242  + diff_search_vector(ispin)%local_data(:, :)
1243  END DO
1244 
1245  ! Make new to old
1246  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  do_restart = (norm_residual >= norm_residual_old .OR. (mod(restart_counter, restart_every) == 0))
1251 
1252  restart_counter = restart_counter + 1
1253  norm_residual_old = norm_residual
1254 
1255  END DO
1256 
1257  IF (unit_nr > 0) THEN
1258  WRITE (unit_nr, '(T4,A)') repeat("-", 40)
1259  IF (converged) THEN
1260  WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations converged in', cycles_passed, ' steps'
1261  ELSE
1262  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  DEALLOCATE (fm_struct_tmp)
1267  CALL cp_fm_release(residual)
1268  CALL cp_fm_release(diff_search_vector)
1269  CALL cp_fm_release(search_vector)
1270  CALL cp_fm_release(a_dot_search_vector)
1271 
1272  CALL timestop(handle)
1273 
1274  END SUBROUTINE solve_z_vector_cg
1275 
1276 ! **************************************************************************************************
1277 !> \brief ...
1278 !> \param matrix1 ...
1279 !> \param matrix2 ...
1280 !> \return ...
1281 ! **************************************************************************************************
1282  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  dotproduct = 0.0_dp
1289  DO ispin = 1, SIZE(matrix1)
1290  dotproduct = dotproduct + accurate_dot_product(matrix1(ispin)%local_data, matrix2(ispin)%local_data)
1291  END DO
1292  CALL matrix1(1)%matrix_struct%para_env%sum(dotproduct)
1293 
1294  END FUNCTION accurate_dot_product_spin
1295 
1296 ! **************************************************************************************************
1297 !> \brief ...
1298 !> \param qs_env ...
1299 ! **************************************************************************************************
1300  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  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  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  TARGET :: matrix_ks_aux
1319  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p_mp2, &
1320  matrix_p_mp2_admm, matrix_s, rho1, &
1321  rho_ao, rho_ao_aux, scrm
1322  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, scrm_kp
1323  TYPE(dft_control_type), POINTER :: dft_control
1324  TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1325  TYPE(linres_control_type), POINTER :: linres_control
1326  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  POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
1330  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1331  TYPE(pw_c1d_gs_type) :: pot_g, rho_tot_g, temp_pw_g
1332  TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: dvg
1333  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  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_mp2_r, rho_mp2_r_aux, rho_r, &
1340  tau_mp2_r, vadmm_rspace, vtau_rspace, &
1341  vxc_rspace
1342  TYPE(qs_dispersion_type), POINTER :: dispersion_env
1343  TYPE(qs_energy_type), POINTER :: energy
1344  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1345  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  CALL timeset(routinen, handle)
1355 
1356  NULLIFY (input, pw_env, matrix_s, rho, energy, force, virial, &
1357  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  x_data=x_data)
1377 
1378  logger => cp_get_default_logger()
1379  iounit = cp_print_key_unit_nr(logger, input, "DFT%XC%WF_CORRELATION%PRINT", &
1380  extension=".mp2Log")
1381 
1382  do_exx = .false.
1383  IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
1384  hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
1385  CALL section_vals_get(hfx_section, explicit=do_exx)
1386  END IF
1387 
1388  nimages = dft_control%nimages
1389  cpassert(nimages == 1)
1390  NULLIFY (cell_to_index)
1391 
1392  p_env => qs_env%mp2_env%ri_grad%p_env
1393 
1394  CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp, rho_r=rho_r, rho_g=rho_g)
1395  nspins = SIZE(rho_ao)
1396 
1397  ! check if we have to calculate the virial
1398  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1399  IF (use_virial) virial%pv_calculate = .true.
1400 
1401  CALL zero_qs_force(force)
1402  IF (use_virial) CALL zero_virial(virial, .false.)
1403 
1404  DO ispin = 1, nspins
1405  CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
1406  END DO
1407 
1408  IF (nspins == 2) THEN
1409  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  NULLIFY (scrm)
1414  IF (debug_forces) THEN
1415  deb(1:3) = force(1)%kinetic(1:3, 1)
1416  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  matrix_p=rho_ao(1)%matrix)
1423  IF (debug_forces) THEN
1424  deb(1:3) = force(1)%kinetic(1:3, 1) - deb(1:3)
1425  CALL para_env%sum(deb)
1426  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", deb
1427  IF (use_virial) THEN
1428  e_dummy = third_tr(virial%pv_virial) - e_dummy
1429  CALL para_env%sum(e_dummy)
1430  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: T ", e_dummy
1431  END IF
1432  END IF
1433  CALL dbcsr_deallocate_matrix_set(scrm)
1434 
1435  IF (nspins == 2) THEN
1436  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  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  atomic_kind_set=atomic_kind_set, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl)
1443  IF (ASSOCIATED(sac_ae)) THEN
1444  IF (debug_forces) THEN
1445  deb(1:3) = force(1)%all_potential(1:3, 1)
1446  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  nimages, cell_to_index)
1452  IF (debug_forces) THEN
1453  deb(1:3) = force(1)%all_potential(1:3, 1) - deb(1:3)
1454  CALL para_env%sum(deb)
1455  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", deb
1456  IF (use_virial) THEN
1457  e_dummy = third_tr(virial%pv_virial) - e_dummy
1458  CALL para_env%sum(e_dummy)
1459  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  IF (ASSOCIATED(sac_ppl)) THEN
1464  IF (debug_forces) THEN
1465  deb(1:3) = force(1)%gth_ppl(1:3, 1)
1466  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  nimages=1, basis_type="ORB")
1472  IF (debug_forces) THEN
1473  deb(1:3) = force(1)%gth_ppl(1:3, 1) - deb(1:3)
1474  CALL para_env%sum(deb)
1475  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", deb
1476  IF (use_virial) THEN
1477  e_dummy = third_tr(virial%pv_virial) - e_dummy
1478  CALL para_env%sum(e_dummy)
1479  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  IF (ASSOCIATED(sap_ppnl)) THEN
1484  IF (debug_forces) THEN
1485  deb(1:3) = force(1)%gth_ppnl(1:3, 1)
1486  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  nimages=1, basis_type="ORB")
1492  IF (debug_forces) THEN
1493  deb(1:3) = force(1)%gth_ppnl(1:3, 1) - deb(1:3)
1494  CALL para_env%sum(deb)
1495  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", deb
1496  IF (use_virial) THEN
1497  e_dummy = third_tr(virial%pv_virial) - e_dummy
1498  CALL para_env%sum(e_dummy)
1499  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  NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)
1506  IF (use_virial) THEN
1507  h_stress = 0.0_dp
1508  CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
1509  ! Update virial
1510  virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1511  virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1512  IF (.NOT. do_exx) THEN
1513  virial%pv_exc = virial%pv_exc - virial%pv_xc
1514  virial%pv_virial = virial%pv_virial - virial%pv_xc
1515  ELSE
1516  virial%pv_xc = 0.0_dp
1517  END IF
1518  IF (debug_forces) THEN
1519  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: xc ", third_tr(h_stress)
1520  CALL para_env%sum(virial%pv_xc(1, 1))
1521  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  CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc)
1525  END IF
1526 
1527  ! Vhxc
1528  CALL get_qs_env(qs_env, pw_env=pw_env)
1529  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1530  poisson_env=poisson_env)
1531  CALL auxbas_pw_pool%create_pw(vhxc_rspace)
1532  IF (use_virial) h_stress = virial%pv_virial
1533  IF (debug_forces) THEN
1534  deb(1:3) = force(1)%rho_elec(1:3, 1)
1535  IF (use_virial) e_dummy = third_tr(h_stress)
1536  END IF
1537  IF (do_exx) THEN
1538  DO ispin = 1, nspins
1539  CALL pw_transfer(vh_rspace, vhxc_rspace)
1540  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  qs_env=qs_env, calculate_forces=.true.)
1544  CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
1545  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  qs_env=qs_env, calculate_forces=.true.)
1549  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  qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
1553  END IF
1554  END DO
1555  ELSE
1556  DO ispin = 1, nspins
1557  CALL pw_transfer(vh_rspace, vhxc_rspace)
1558  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  qs_env=qs_env, calculate_forces=.true.)
1562  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  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  deb(1:3) = force(1)%rho_elec(1:3, 1) - deb(1:3)
1571  CALL para_env%sum(deb)
1572  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dVhxc ", deb
1573  IF (use_virial) THEN
1574  e_dummy = third_tr(virial%pv_virial) - e_dummy
1575  CALL para_env%sum(e_dummy)
1576  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Vhxc ", e_dummy
1577  END IF
1578  END IF
1579  IF (use_virial) THEN
1580  h_stress = virial%pv_virial - h_stress
1581  virial%pv_ehartree = virial%pv_ehartree + h_stress
1582 
1583  CALL qs_rho_get(p_env%rho1, rho_r=rho_mp2_r, tau_r=tau_mp2_r)
1584  e_xc = 0.0_dp
1585  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  CALL pw_scale(vxc_rspace(ispin), 1.0_dp/vxc_rspace(ispin)%pw_grid%dvol)
1588  e_xc = e_xc + pw_integral_ab(rho_mp2_r(ispin), vxc_rspace(ispin))
1589  IF (ASSOCIATED(vtau_rspace)) CALL pw_scale(vtau_rspace(ispin), 1.0_dp/vtau_rspace(ispin)%pw_grid%dvol)
1590  IF (ASSOCIATED(vtau_rspace)) e_xc = e_xc + pw_integral_ab(tau_mp2_r(ispin), vtau_rspace(ispin))
1591  END DO
1592  IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG VIRIAL:: vxc*d1 ", e_xc
1593  DO alpha = 1, 3
1594  virial%pv_exc(alpha, alpha) = virial%pv_exc(alpha, alpha) - e_xc/real(para_env%num_pe, dp)
1595  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  DO ispin = 1, nspins
1599  CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
1600  IF (ASSOCIATED(vtau_rspace)) THEN
1601  CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
1602  END IF
1603  END DO
1604  DEALLOCATE (vxc_rspace)
1605  CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
1606  IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
1607 
1608  DO ispin = 1, nspins
1609  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  NULLIFY (poisson_env, auxbas_pw_pool)
1614  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1615  poisson_env=poisson_env)
1616 
1617  ! get some of the grids ready
1618  CALL auxbas_pw_pool%create_pw(pot_r)
1619  CALL auxbas_pw_pool%create_pw(pot_g)
1620  CALL auxbas_pw_pool%create_pw(rho_tot_g)
1621 
1622  CALL pw_zero(rho_tot_g)
1623 
1624  CALL qs_rho_get(p_env%rho1, rho_r=rho_mp2_r, rho_g=rho_mp2_g)
1625  DO ispin = 1, nspins
1626  CALL pw_axpy(rho_mp2_g(ispin), rho_tot_g)
1627  END DO
1628 
1629  IF (use_virial) THEN
1630  ALLOCATE (dvg(3))
1631  DO idir = 1, 3
1632  CALL auxbas_pw_pool%create_pw(dvg(idir))
1633  END DO
1634  CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g, dvhartree=dvg)
1635  ELSE
1636  CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g)
1637  END IF
1638 
1639  CALL pw_transfer(pot_g, pot_r)
1640  CALL pw_scale(pot_r, pot_r%pw_grid%dvol)
1641  CALL pw_axpy(pot_r, vh_rspace)
1642 
1643  ! calculate core forces
1644  CALL integrate_v_core_rspace(vh_rspace, qs_env)
1645  IF (debug_forces) THEN
1646  deb(:) = force(1)%rho_core(:, 1)
1647  CALL para_env%sum(deb)
1648  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: core ", deb
1649  IF (use_virial) THEN
1650  e_dummy = third_tr(virial%pv_virial) - e_dummy
1651  CALL para_env%sum(e_dummy)
1652  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: core ", e_dummy
1653  END IF
1654  END IF
1655  CALL auxbas_pw_pool%give_back_pw(vh_rspace)
1656 
1657  IF (use_virial) THEN
1658  ! update virial if necessary with the volume term
1659  ! first create pw auxiliary stuff
1660  CALL auxbas_pw_pool%create_pw(temp_pw_g)
1661 
1662  ! make a copy of the MP2 density in G space
1663  CALL pw_copy(rho_tot_g, temp_pw_g)
1664 
1665  ! calculate total SCF density and potential
1666  CALL pw_copy(rho_g(1), rho_tot_g)
1667  IF (nspins == 2) CALL pw_axpy(rho_g(2), rho_tot_g)
1668  CALL pw_axpy(rho_core, rho_tot_g)
1669  CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g)
1670 
1671  ! finally update virial with the volume contribution
1672  e_hartree = pw_integral_ab(temp_pw_g, pot_g)
1673  IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: vh1*d0 ", e_hartree
1674 
1675  h_stress = 0.0_dp
1676  DO alpha = 1, 3
1677  comp = 0
1678  comp(alpha) = 1
1679  CALL pw_copy(pot_g, rho_tot_g)
1680  CALL pw_derive(rho_tot_g, comp)
1681  h_stress(alpha, alpha) = -e_hartree
1682  DO beta = alpha, 3
1683  h_stress(alpha, beta) = h_stress(alpha, beta) &
1684  - 2.0_dp*pw_integral_ab(rho_tot_g, dvg(beta))/fourpi
1685  h_stress(beta, alpha) = h_stress(alpha, beta)
1686  END DO
1687  END DO
1688  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  CALL auxbas_pw_pool%give_back_pw(temp_pw_g)
1692  DO idir = 1, 3
1693  CALL auxbas_pw_pool%give_back_pw(dvg(idir))
1694  END DO
1695  DEALLOCATE (dvg)
1696 
1697  virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1698  virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1699 
1700  END IF
1701 
1702  DO ispin = 1, nspins
1703  CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
1704  IF (dft_control%do_admm) CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
1705  END DO
1706 
1707  CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
1708 
1709  IF (dft_control%do_admm) THEN
1710  CALL get_qs_env(qs_env, admm_env=admm_env)
1711  xc_section => admm_env%xc_section_primary
1712  ELSE
1713  xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1714  END IF
1715 
1716  IF (use_virial) THEN
1717  h_stress = 0.0_dp
1718  pv_virial = virial%pv_virial
1719  END IF
1720  IF (debug_forces) THEN
1721  deb = force(1)%rho_elec(1:3, 1)
1722  IF (use_virial) e_dummy = third_tr(pv_virial)
1723  END IF
1724  CALL apply_2nd_order_kernel(qs_env, p_env, .false., .true., use_virial, h_stress)
1725  IF (use_virial) THEN
1726  virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_virial)
1727  IF (debug_forces) THEN
1728  e_dummy = third_tr(virial%pv_virial - pv_virial)
1729  CALL para_env%sum(e_dummy)
1730  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Kh ", e_dummy
1731  END IF
1732  virial%pv_exc = virial%pv_exc + h_stress
1733  virial%pv_virial = virial%pv_virial + h_stress
1734  IF (debug_forces) THEN
1735  e_dummy = third_tr(h_stress)
1736  CALL para_env%sum(e_dummy)
1737  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  deb(:) = force(1)%rho_elec(:, 1) - deb
1742  CALL para_env%sum(deb)
1743  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P0*Khxc ", deb
1744  IF (use_virial) THEN
1745  e_dummy = third_tr(virial%pv_virial) - e_dummy
1746  CALL para_env%sum(e_dummy)
1747  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  NULLIFY (hfx_sections)
1753  hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1754  CALL section_vals_get(hfx_sections, explicit=do_hfx)
1755  IF (do_hfx) THEN
1756  IF (do_exx) THEN
1757  IF (dft_control%do_admm) THEN
1758  CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
1759  CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux, rho_ao_kp=rho_ao_kp)
1760  rho1 => p_env%p1_admm
1761  ELSE
1762  rho1 => p_env%p1
1763  END IF
1764  ELSE
1765  IF (dft_control%do_admm) THEN
1766  CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
1767  CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux, rho_ao_kp=rho_ao_kp)
1768  DO ispin = 1, nspins
1769  CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
1770  END DO
1771  rho1 => p_env%p1_admm
1772  ELSE
1773  DO ispin = 1, nspins
1774  CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
1775  END DO
1776  rho1 => p_env%p1
1777  END IF
1778  END IF
1779 
1780  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  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  1, use_virial, resp_only=do_exx)
1790  END IF
1791 
1792  IF (use_virial) THEN
1793  virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1794  virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1795  END IF
1796  IF (debug_forces) THEN
1797  deb(1:3) = force(1)%fock_4c(1:3, 1) - deb(1:3)
1798  CALL para_env%sum(deb)
1799  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx ", deb
1800  IF (use_virial) THEN
1801  e_dummy = third_tr(virial%pv_fock_4c)
1802  CALL para_env%sum(e_dummy)
1803  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: hfx ", e_dummy
1804  END IF
1805  END IF
1806 
1807  IF (.NOT. do_exx) THEN
1808  IF (dft_control%do_admm) THEN
1809  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1810  DO ispin = 1, nspins
1811  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  DO ispin = 1, nspins
1815  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  IF (dft_control%do_admm) THEN
1821  IF (debug_forces) THEN
1822  deb = force(1)%overlap_admm(1:3, 1)
1823  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  IF (nspins == 1) CALL dbcsr_scale(p_env%kpp1_admm(1)%matrix, 0.5_dp)
1827  CALL admm_projection_derivative(qs_env, p_env%kpp1_admm, rho_ao)
1828  IF (debug_forces) THEN
1829  deb(:) = force(1)%overlap_admm(:, 1) - deb
1830  CALL para_env%sum(deb)
1831  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*KADMM*dS'", deb
1832  IF (use_virial) THEN
1833  e_dummy = third_tr(virial%pv_virial) - e_dummy
1834  CALL para_env%sum(e_dummy)
1835  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: KADMM*S' ", e_dummy
1836  END IF
1837  END IF
1838 
1839  ALLOCATE (matrix_ks_aux(nspins))
1840  DO ispin = 1, nspins
1841  NULLIFY (matrix_ks_aux(ispin)%matrix)
1842  ALLOCATE (matrix_ks_aux(ispin)%matrix)
1843  CALL dbcsr_copy(matrix_ks_aux(ispin)%matrix, p_env%kpp1_admm(ispin)%matrix)
1844  CALL dbcsr_set(matrix_ks_aux(ispin)%matrix, 0.0_dp)
1845  END DO
1846 
1847  ! Calculate kernel
1848  CALL tddft_hfx_matrix(matrix_ks_aux, rho_ao_aux, qs_env, .false.)
1849 
1850  IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
1851  CALL get_qs_env(qs_env, ks_env=ks_env)
1852  CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list_aux_fit)
1853 
1854  DO ispin = 1, nspins
1855  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  IF (use_virial) THEN
1859  CALL qs_rho_get(p_env%rho1_admm, rho_r=rho_mp2_r_aux)
1860  e_xc = 0.0_dp
1861  DO ispin = 1, nspins
1862  e_xc = e_xc + pw_integral_ab(rho_mp2_r_aux(ispin), vadmm_rspace(ispin))
1863  END DO
1864 
1865  e_xc = -e_xc/vadmm_rspace(1)%pw_grid%dvol/real(para_env%num_pe, dp)
1866 
1867  ! Update the virial
1868  DO alpha = 1, 3
1869  virial%pv_exc(alpha, alpha) = virial%pv_exc(alpha, alpha) + e_xc
1870  virial%pv_virial(alpha, alpha) = virial%pv_virial(alpha, alpha) + e_xc
1871  END DO
1872  IF (debug_forces) THEN
1873  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: P1*VADMM ", e_xc
1874  END IF
1875  END IF
1876 
1877  IF (use_virial) h_stress = virial%pv_virial
1878  IF (debug_forces) THEN
1879  deb = force(1)%rho_elec(1:3, 1)
1880  IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1881  END IF
1882  DO ispin = 1, nspins
1883  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  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  task_list_external=task_list_aux_fit, pmat=rho_ao_aux(ispin))
1891  END IF
1892  CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
1893  END DO
1894  IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - h_stress)
1895  DEALLOCATE (vadmm_rspace)
1896  IF (debug_forces) THEN
1897  deb(:) = force(1)%rho_elec(:, 1) - deb
1898  CALL para_env%sum(deb)
1899  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*VADMM' ", deb
1900  IF (use_virial) THEN
1901  e_dummy = third_tr(virial%pv_virial) - e_dummy
1902  CALL para_env%sum(e_dummy)
1903  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: VADMM' ", e_dummy
1904  END IF
1905  END IF
1906 
1907  DO ispin = 1, nspins
1908  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  DO ispin = 1, nspins
1914  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  deb = force(1)%overlap_admm(1:3, 1)
1919  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  IF (do_exx) THEN
1923  CALL admm_projection_derivative(qs_env, matrix_ks_aux, matrix_p_mp2)
1924  ELSE
1925  CALL admm_projection_derivative(qs_env, matrix_ks_aux, rho_ao)
1926  END IF
1927  IF (debug_forces) THEN
1928  deb(:) = force(1)%overlap_admm(:, 1) - deb
1929  CALL para_env%sum(deb)
1930  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*VADMM*dS'", deb
1931  IF (use_virial) THEN
1932  e_dummy = third_tr(virial%pv_virial) - e_dummy
1933  CALL para_env%sum(e_dummy)
1934  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: VADMM*S' ", e_dummy
1935  END IF
1936  END IF
1937 
1938  DO ispin = 1, nspins
1939  CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, -1.0_dp)
1940  END DO
1941 
1942  DO ispin = 1, nspins
1943  CALL dbcsr_release(matrix_ks_aux(ispin)%matrix)
1944  DEALLOCATE (matrix_ks_aux(ispin)%matrix)
1945  END DO
1946  DEALLOCATE (matrix_ks_aux)
1947  END IF
1948  END IF
1949 
1950  CALL dbcsr_scale(p_env%w1(1)%matrix, -1.0_dp)
1951 
1952  ! Finish matrix_w_mp2 with occ-occ block
1953  DO ispin = 1, nspins
1954  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  p_env%w1(1)%matrix, 1.0_dp, nocc)
1957  END DO
1958 
1959  IF (debug_forces .AND. use_virial) e_dummy = third_tr(virial%pv_virial)
1960 
1961  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  matrix_p=p_env%w1(1)%matrix)
1967  CALL dbcsr_deallocate_matrix_set(scrm)
1968 
1969  IF (debug_forces) THEN
1970  deb = force(1)%overlap(1:3, 1)
1971  CALL para_env%sum(deb)
1972  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: W*dS ", deb
1973  IF (use_virial) THEN
1974  e_dummy = third_tr(virial%pv_virial) - e_dummy
1975  CALL para_env%sum(e_dummy)
1976  IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: S ", e_dummy
1977  END IF
1978  END IF
1979 
1980  CALL auxbas_pw_pool%give_back_pw(pot_r)
1981  CALL auxbas_pw_pool%give_back_pw(pot_g)
1982  CALL auxbas_pw_pool%give_back_pw(rho_tot_g)
1983 
1984  ! Release linres stuff
1985  CALL p_env_release(p_env)
1986  DEALLOCATE (p_env)
1987  NULLIFY (qs_env%mp2_env%ri_grad%p_env)
1988 
1989  CALL sum_qs_force(force, qs_env%mp2_env%ri_grad%mp2_force)
1990  CALL deallocate_qs_force(qs_env%mp2_env%ri_grad%mp2_force)
1991 
1992  IF (use_virial) THEN
1993  virial%pv_mp2 = qs_env%mp2_env%ri_grad%mp2_virial
1994  virial%pv_virial = virial%pv_virial + qs_env%mp2_env%ri_grad%mp2_virial
1995  IF (debug_forces) THEN
1996  e_dummy = third_tr(virial%pv_mp2)
1997  CALL para_env%sum(e_dummy)
1998  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  IF (use_virial) virial%pv_calculate = .false.
2003 
2004  ! Add the dispersion forces and virials
2005  CALL get_qs_env(qs_env, dispersion_env=dispersion_env)
2006  CALL calculate_dispersion_pairpot(qs_env, dispersion_env, e_dummy, .true.)
2007 
2008  CALL cp_print_key_finished_output(iounit, logger, input, &
2009  "DFT%XC%WF_CORRELATION%PRINT")
2010 
2011  CALL timestop(handle)
2012 
2013  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  PURE FUNCTION third_tr(matrix)
2021  REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: matrix
2022  REAL(kind=dp) :: third_tr
2023 
2024  third_tr = (matrix(1, 1) + matrix(2, 2) + matrix(3, 3))/3.0_dp
2025 
2026  END FUNCTION third_tr
2027 
2028 END MODULE mp2_cphf
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
admm_methods
Contains ADMM methods which require molecular orbitals.
Definition: admm_methods.F:15
admm_methods::admm_projection_derivative
subroutine, public admm_projection_derivative(qs_env, matrix_hz, matrix_pz, fval)
Calculate derivatives terms from overlap matrices.
Definition: admm_methods.F:3239
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition: admm_types.F:593
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
core_ae
Calculation of the nuclear attraction contribution to the core Hamiltonian <a|erfc|b> :we only calcul...
Definition: core_ae.F:14
core_ae::build_core_ae
subroutine, public build_core_ae(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, nimages, cell_to_index)
...
Definition: core_ae.F:84
core_ppl
Calculation of the local pseudopotential contribution to the core Hamiltonian <a|V(local)|b> = <a|Sum...
Definition: core_ppl.F:18
core_ppl::build_core_ppl
subroutine, public build_core_ppl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, nimages, cell_to_index, basis_type, deltaR)
...
Definition: core_ppl.F:95
core_ppnl
Calculation of the non-local pseudopotential contribution to the core Hamiltonian <a|V(non-local)|b> ...
Definition: core_ppnl.F:15
core_ppnl::build_core_ppnl
subroutine, public build_core_ppnl(matrix_h, matrix_p, force, virial, calculate_forces, use_virial, nder, qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, eps_ppnl, nimages, cell_to_index, basis_type, deltaR, matrix_l)
...
Definition: core_ppnl.F:89
cp_blacs_env
methods related to the blacs parallel environment
Definition: cp_blacs_env.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition: cp_dbcsr_operations.F:887
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition: cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_upper_to_full
subroutine, public cp_fm_upper_to_full(matrix, work)
given an upper triangular matrix computes the corresponding full matrix
Definition: cp_fm_basic_linalg.F:1693
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_to_fm_submat
subroutine, public cp_fm_to_fm_submat(msource, mtarget, nrow, ncol, s_firstrow, s_firstcol, t_firstrow, t_firstcol)
copy just a part ot the matrix
Definition: cp_fm_types.F:1473
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
hfx_admm_utils
Utilities for hfx and admm methods.
Definition: hfx_admm_utils.F:17
hfx_admm_utils::tddft_hfx_matrix
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data)
Add the hfx contributions to the Hamiltonian.
Definition: hfx_admm_utils.F:2477
hfx_derivatives
Routines to calculate derivatives with respect to basis function origin.
Definition: hfx_derivatives.F:14
hfx_derivatives::derivatives_four_center
subroutine, public derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, irep, use_virial, adiabatic_rescale_factor, resp_only, external_x_data)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Definition: hfx_derivatives.F:117
hfx_exx
Routines to calculate EXX in RPA and energy correction methods.
Definition: hfx_exx.F:16
hfx_exx::add_exx_to_rhs
subroutine, public add_exx_to_rhs(rhs, qs_env, ext_hfx_section, x_data, recalc_integrals, do_admm, do_ec, do_exx, reuse_hfx)
Add the EXX contribution to the RHS of the Z-vector equation, namely the HF Hamiltonian.
Definition: hfx_exx.F:325
hfx_ri
RI-methods for HFX.
Definition: hfx_ri.F:12
hfx_ri::hfx_ri_update_forces
subroutine, public hfx_ri_update_forces(qs_env, ri_data, nspins, hf_fraction, rho_ao, rho_ao_resp, mos, use_virial, resp_only, rescale_factor)
the general routine that calls the relevant force code
Definition: hfx_ri.F:3036
hfx_types
Types and set/get functions for HFX.
Definition: hfx_types.F:15
hfx_types::alloc_containers
subroutine, public alloc_containers(DATA, bin_size)
...
Definition: hfx_types.F:2906
hfx_types::hfx_init_container
subroutine, public hfx_init_container(container, memory_usage, do_disk_storage)
This routine deletes all list entries in a container in order to deallocate the memory.
Definition: hfx_types.F:2523
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::z_solver_pople
integer, parameter, public z_solver_pople
Definition: input_constants.F:1057
input_constants::z_solver_cg
integer, parameter, public z_solver_cg
Definition: input_constants.F:1057
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition: input_constants.F:863
input_constants::z_solver_sd
integer, parameter, public z_solver_sd
Definition: input_constants.F:1057
input_constants::z_solver_richardson
integer, parameter, public z_solver_richardson
Definition: input_constants.F:1057
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition: input_constants.F:513
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition: input_section_types.F:697
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition: kahan_sum.F:29
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
linear_systems
Provides interfaces to LAPACK routines for factorisation and linear system solving.
Definition: linear_systems.F:18
linear_systems::solve_system
subroutine, public solve_system(matrix, mysize, eigenvectors)
...
Definition: linear_systems.F:41
machine
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
machine::m_walltime
real(kind=dp) function, public m_walltime()
returns time from a real-time clock, protected against rolling early/easily
Definition: machine.F:123
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition: mathconstants.F:113
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
mp2_cphf
Routines to calculate CPHF like update and solve Z-vector equation for MP2 gradients (only GPW)
Definition: mp2_cphf.F:14
mp2_cphf::solve_z_vector_eq
subroutine, public solve_z_vector_eq(qs_env, mp2_env, para_env, dft_control, mo_coeff, nmo, homo, Eigenval, unit_nr)
Solve Z-vector equations necessary for the calculation of the MP2 gradients, in order to be consisten...
Definition: mp2_cphf.F:161
mp2_cphf::update_mp2_forces
subroutine, public update_mp2_forces(qs_env)
...
Definition: mp2_cphf.F:1301
mp2_types
Types needed for MP2 calculations.
Definition: mp2_types.F:14
parallel_gemm_api
basic linear algebra operations for full matrixes
Definition: parallel_gemm_api.F:14
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_methods::pw_derive
subroutine, public pw_derive(pw, n)
Calculate the derivative of a plane wave vector.
Definition: pw_methods.F:10106
pw_poisson_methods
Definition: pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_2nd_kernel_ao
Routines to calculate 2nd order kernels from a given response density in ao basis linear response scf...
Definition: qs_2nd_kernel_ao.F:15
qs_2nd_kernel_ao::apply_2nd_order_kernel
subroutine, public apply_2nd_order_kernel(qs_env, p_env, recalc_hfx_integrals, calc_forces, calc_virial, virial)
Calculate a second order kernel (DFT, HF, ADMM correction) for a given density.
Definition: qs_2nd_kernel_ao.F:113
qs_density_matrices
collects routines that calculate density matrices
Definition: qs_density_matrices.F:14
qs_density_matrices::calculate_whz_matrix
subroutine, public calculate_whz_matrix(c0vec, hzm, w_matrix, focc, nocc)
Calculate the Wz matrix from the MO eigenvectors, MO eigenvalues, and the MO occupation numbers....
Definition: qs_density_matrices.F:425
qs_dispersion_pairpot
Calculation of dispersion using pair potentials.
Definition: qs_dispersion_pairpot.F:12
qs_dispersion_pairpot::calculate_dispersion_pairpot
subroutine, public calculate_dispersion_pairpot(qs_env, dispersion_env, energy, calculate_forces, atevdw)
...
Definition: qs_dispersion_pairpot.F:1245
qs_dispersion_types
Definition of disperson types for DFT calculations.
Definition: qs_dispersion_types.F:12
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, WannierCentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Definition: qs_environment_types.F:1035
qs_force_types
Definition: qs_force_types.F:13
qs_force_types::sum_qs_force
subroutine, public sum_qs_force(qs_force_out, qs_force_in)
Sum up two qs_force entities qs_force_out = qs_force_out + qs_force_in.
Definition: qs_force_types.F:292
qs_force_types::deallocate_qs_force
subroutine, public deallocate_qs_force(qs_force)
Deallocate a Quickstep force data structure.
Definition: qs_force_types.F:131
qs_force_types::zero_qs_force
subroutine, public zero_qs_force(qs_force)
Initialize a Quickstep force data structure.
Definition: qs_force_types.F:250
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition: qs_integrate_potential.F:22
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kinetic
Calculation of kinetic energy matrix and forces.
Definition: qs_kinetic.F:15
qs_kinetic::build_kinetic_matrix
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
Definition: qs_kinetic.F:101
qs_ks_reference
Calculate the KS reference potentials.
Definition: qs_ks_reference.F:14
qs_ks_reference::ks_ref_potential
subroutine, public ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
calculate the Kohn-Sham reference potential
Definition: qs_ks_reference.F:95
qs_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition: qs_ks_types.F:567
qs_linres_types
Type definitiona for linear response calculations.
Definition: qs_linres_types.F:12
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_overlap
Calculation of overlap matrix, its derivatives and forces.
Definition: qs_overlap.F:19
qs_overlap::build_overlap_matrix
subroutine, public build_overlap_matrix(ks_env, matrix_s, matrixkp_s, matrix_name, nderivative, basis_type_a, basis_type_b, sab_nl, calculate_forces, matrix_p, matrixkp_p)
Calculation of the overlap matrix over Cartesian Gaussian functions.
Definition: qs_overlap.F:120
qs_p_env_methods
Utility functions for the perturbation calculations.
Definition: qs_p_env_methods.F:16
qs_p_env_methods::p_env_psi0_changed
subroutine, public p_env_psi0_changed(p_env, qs_env)
To be called after the value of psi0 has changed. Recalculates the quantities S_psi0 and m_epsilon.
Definition: qs_p_env_methods.F:482
qs_p_env_methods::p_env_create
subroutine, public p_env_create(p_env, qs_env, p1_option, p1_admm_option, orthogonal_orbitals, linres_control)
allocates and initializes the perturbation environment (no setup)
Definition: qs_p_env_methods.F:133
qs_p_env_methods::p_env_update_rho
subroutine, public p_env_update_rho(p_env, qs_env)
...
Definition: qs_p_env_methods.F:398
qs_p_env_methods::p_env_check_i_alloc
subroutine, public p_env_check_i_alloc(p_env, qs_env)
checks that the intenal storage is allocated, and allocs it if needed
Definition: qs_p_env_methods.F:317
qs_p_env_types
basis types for the calculation of the perturbation of density theory.
Definition: qs_p_env_types.F:14
qs_p_env_types::p_env_release
subroutine, public p_env_release(p_env)
relases the given p_env (see doc/ReferenceCounting.html)
Definition: qs_p_env_types.F:99
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
task_list_types
types for task lists
Definition: task_list_types.F:14
virial_types
Definition: virial_types.F:13
virial_types::zero_virial
subroutine, public zero_virial(virial, reset)
...
Definition: virial_types.F:123