(git:ed6f26b)
Loading...
Searching...
No Matches
mp2_cphf.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Routines to calculate CPHF like update and solve Z-vector equation
10!> for MP2 gradients (only GPW)
11!> \par History
12!> 11.2013 created [Mauro Del Ben]
13! **************************************************************************************************
14MODULE mp2_cphf
15 USE admm_methods, ONLY: admm_projection_derivative
16 USE admm_types, ONLY: admm_type,&
17 get_admm_env
18 USE atomic_kind_types, ONLY: atomic_kind_type
19 USE cell_types, ONLY: cell_type
20 USE core_ae, ONLY: build_core_ae
21 USE core_ppl, ONLY: build_core_ppl
22 USE core_ppnl, ONLY: build_core_ppnl
23 USE cp_blacs_env, ONLY: cp_blacs_env_type
24 USE cp_control_types, ONLY: dft_control_type
25 USE cp_dbcsr_api, ONLY: dbcsr_add,&
26 dbcsr_copy,&
27 dbcsr_p_type,&
28 dbcsr_release,&
29 dbcsr_scale,&
30 dbcsr_set
31 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
32 cp_dbcsr_plus_fm_fm_t,&
33 dbcsr_allocate_matrix_set,&
34 dbcsr_deallocate_matrix_set
35 USE cp_fm_basic_linalg, ONLY: cp_fm_uplo_to_full
36 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37 cp_fm_struct_p_type,&
38 cp_fm_struct_release,&
39 cp_fm_struct_type
40 USE cp_fm_types, ONLY: cp_fm_create,&
41 cp_fm_get_info,&
42 cp_fm_release,&
43 cp_fm_set_all,&
44 cp_fm_to_fm_submat,&
45 cp_fm_type
46 USE cp_log_handling, ONLY: cp_get_default_logger,&
47 cp_logger_type
48 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
49 cp_print_key_unit_nr
50 USE hfx_admm_utils, ONLY: tddft_hfx_matrix
51 USE hfx_derivatives, ONLY: derivatives_four_center
52 USE hfx_exx, ONLY: add_exx_to_rhs
53 USE hfx_ri, ONLY: hfx_ri_update_forces
54 USE hfx_types, ONLY: alloc_containers,&
55 hfx_container_type,&
56 hfx_init_container,&
57 hfx_type
58 USE input_constants, ONLY: do_admm_aux_exch_func_none,&
59 ot_precond_full_all,&
60 z_solver_cg,&
61 z_solver_pople,&
62 z_solver_richardson,&
63 z_solver_sd
64 USE input_section_types, ONLY: section_vals_get,&
65 section_vals_get_subs_vals,&
66 section_vals_type
67 USE kahan_sum, ONLY: accurate_dot_product
68 USE kinds, ONLY: dp
69 USE linear_systems, ONLY: solve_system
70 USE machine, ONLY: m_flush,&
71 m_walltime
72 USE mathconstants, ONLY: fourpi
73 USE message_passing, ONLY: mp_para_env_type
74 USE mp2_types, ONLY: mp2_type,&
75 ri_rpa_method_gpw
76 USE parallel_gemm_api, ONLY: parallel_gemm
77 USE particle_types, ONLY: particle_type
78 USE pw_env_types, ONLY: pw_env_get,&
79 pw_env_type
80 USE pw_methods, ONLY: pw_axpy,&
81 pw_copy,&
82 pw_derive,&
83 pw_integral_ab,&
84 pw_scale,&
85 pw_transfer,&
86 pw_zero
87 USE pw_poisson_methods, ONLY: pw_poisson_solve
88 USE pw_poisson_types, ONLY: pw_poisson_type
89 USE pw_pool_types, ONLY: pw_pool_type
90 USE pw_types, ONLY: pw_c1d_gs_type,&
91 pw_r3d_rs_type
92 USE qs_2nd_kernel_ao, ONLY: apply_2nd_order_kernel
93 USE qs_density_matrices, ONLY: calculate_whz_matrix
94 USE qs_dispersion_pairpot, ONLY: calculate_dispersion_pairpot
95 USE qs_dispersion_types, ONLY: qs_dispersion_type
96 USE qs_energy_types, ONLY: qs_energy_type
97 USE qs_environment_types, ONLY: get_qs_env,&
98 qs_environment_type,&
99 set_qs_env
100 USE qs_force_types, ONLY: deallocate_qs_force,&
101 qs_force_type,&
102 sum_qs_force,&
103 zero_qs_force
104 USE qs_integrate_potential, ONLY: integrate_v_core_rspace,&
105 integrate_v_rspace
106 USE qs_kind_types, ONLY: qs_kind_type
107 USE qs_kinetic, ONLY: build_kinetic_matrix
108 USE qs_ks_reference, ONLY: ks_ref_potential
109 USE qs_ks_types, ONLY: qs_ks_env_type,&
110 set_ks_env
111 USE qs_linres_types, ONLY: linres_control_type
112 USE qs_mo_types, ONLY: get_mo_set,&
113 mo_set_type
114 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
115 USE qs_overlap, ONLY: build_overlap_matrix
116 USE qs_p_env_methods, ONLY: p_env_check_i_alloc,&
117 p_env_create,&
118 p_env_psi0_changed,&
119 p_env_update_rho
120 USE qs_p_env_types, ONLY: p_env_release,&
121 qs_p_env_type
122 USE qs_rho_types, ONLY: qs_rho_get,&
123 qs_rho_type
124 USE task_list_types, ONLY: task_list_type
125 USE virial_types, ONLY: virial_type,&
126 zero_virial
127
128!$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
129
130#include "./base/base_uses.f90"
131
132 IMPLICIT NONE
133
134 PRIVATE
135
136 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mp2_cphf'
137 LOGICAL, PARAMETER, PRIVATE :: debug_forces = .true.
138
139 PUBLIC :: solve_z_vector_eq, update_mp2_forces
140
141CONTAINS
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_uplo_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(INOUT) :: fm_ao
561 TYPE(cp_fm_type), INTENT(IN) :: fm_back
562 INTEGER, INTENT(IN) :: transf_type_in
563 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: fm_mo_out
564 LOGICAL, INTENT(IN), OPTIONAL :: recalc_hfx_integrals
565
566 CHARACTER(LEN=*), PARAMETER :: routinen = 'cphf_like_update'
567
568 INTEGER :: handle, i_global, iib, ispin, j_global, &
569 jjb, ncol_local, nrow_local
570 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
571 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho1_ao
572
573 CALL timeset(routinen, handle)
574
575 CALL qs_rho_get(p_env%rho1, rho_ao=rho1_ao)
576
577 ! Determine the first-order density matrices in AO basis
578 DO ispin = 1, nspins
579 CALL dbcsr_set(p_env%p1(ispin)%matrix, 0.0_dp)
580
581 associate(mat_in => fm_mo(ispin))
582
583 ! perform back transformation
584 SELECT CASE (transf_type_in)
585 CASE (1)
586 ! occ-occ block
587 CALL parallel_gemm('N', 'N', dimen, homo(ispin), homo(ispin), 1.0_dp, &
588 mo_coeff_o(ispin), mat_in, 0.0_dp, fm_back)
589 CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_o(ispin), homo(ispin), 1.0_dp, .true.)
590 ! virt-virt block
591 CALL parallel_gemm('N', 'N', dimen, virtual(ispin), virtual(ispin), 1.0_dp, &
592 mo_coeff_v(ispin), mat_in, 0.0_dp, fm_back, &
593 b_first_col=homo(ispin) + 1, &
594 b_first_row=homo(ispin) + 1)
595 CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_v(ispin), virtual(ispin), 1.0_dp, .true.)
596
597 CASE (3)
598 ! virt-occ blocks
599 CALL parallel_gemm('N', 'N', dimen, virtual(ispin), homo(ispin), 1.0_dp, &
600 mo_coeff_o(ispin), mat_in, 0.0_dp, fm_back)
601 CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_v(ispin), virtual(ispin), 1.0_dp, .true.)
602 ! and symmetrize (here again multiply instead of transposing)
603 CALL parallel_gemm('N', 'T', dimen, homo(ispin), virtual(ispin), 1.0_dp, &
604 mo_coeff_v(ispin), mat_in, 0.0_dp, fm_back)
605 CALL cp_dbcsr_plus_fm_fm_t(p_env%p1(ispin)%matrix, fm_back, mo_coeff_o(ispin), homo(ispin), 1.0_dp, .true.)
606
607 CASE DEFAULT
608 ! nothing
609 END SELECT
610 END associate
611
612 CALL dbcsr_copy(rho1_ao(ispin)%matrix, p_env%p1(ispin)%matrix)
613 END DO
614
615 CALL p_env_update_rho(p_env, qs_env)
616
617 IF (transf_type_in == 3) THEN
618 DO ispin = 1, nspins
619
620 ! scale for the orbital energy differences for the diagonal elements
621 CALL cp_fm_get_info(matrix=fm_mo_out(ispin), &
622 nrow_local=nrow_local, &
623 ncol_local=ncol_local, &
624 row_indices=row_indices, &
625 col_indices=col_indices)
626 DO jjb = 1, ncol_local
627 j_global = col_indices(jjb)
628 DO iib = 1, nrow_local
629 i_global = row_indices(iib)
630 fm_mo_out(ispin)%local_data(iib, jjb) = fm_mo(ispin)%local_data(iib, jjb)* &
631 (eigenval(j_global + homo(ispin), ispin) - eigenval(i_global, ispin))
632 END DO
633 END DO
634 END DO
635 END IF
636
637 CALL apply_2nd_order_kernel(qs_env, p_env, recalc_hfx_integrals)
638
639 DO ispin = 1, nspins
640 ! copy back to fm
641 CALL cp_fm_set_all(fm_ao, 0.0_dp)
642 CALL copy_dbcsr_to_fm(matrix=p_env%kpp1(ispin)%matrix, fm=fm_ao)
643 CALL cp_fm_set_all(fm_back, 0.0_dp)
644 CALL cp_fm_uplo_to_full(fm_ao, fm_back)
645
646 associate(mat_out => fm_mo_out(ispin))
647
648 ! transform to MO basis, here we always sum the result into the input matrix
649
650 ! occ-virt block
651 CALL parallel_gemm('T', 'N', homo(ispin), dimen, dimen, 1.0_dp, &
652 mo_coeff_o(ispin), fm_ao, 0.0_dp, fm_back)
653 CALL parallel_gemm('N', 'N', homo(ispin), virtual(ispin), dimen, 1.0_dp, &
654 fm_back, mo_coeff_v(ispin), 1.0_dp, mat_out)
655 END associate
656 END DO
657
658 CALL timestop(handle)
659
660 END SUBROUTINE cphf_like_update
661
662! **************************************************************************************************
663!> \brief Low level subroutine for the iterative solution of a large
664!> system of linear equation
665!> \param qs_env ...
666!> \param mp2_env ...
667!> \param homo ...
668!> \param virtual ...
669!> \param dimen ...
670!> \param unit_nr ...
671!> \param nspins ...
672!> \param mo_coeff_o ...
673!> \param mo_coeff_v ...
674!> \param Eigenval ...
675!> \param p_env ...
676!> \param L_jb ...
677!> \param fm_G_mu_nu ...
678!> \param fm_back ...
679!> \param P_ia ...
680!> \author Mauro Del Ben, Vladimir Rybkin
681! **************************************************************************************************
682 SUBROUTINE solve_z_vector_eq_low(qs_env, mp2_env, homo, virtual, dimen, unit_nr, nspins, &
683 mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
684 L_jb, fm_G_mu_nu, fm_back, P_ia)
685 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
686 TYPE(mp2_type), INTENT(IN) :: mp2_env
687 INTEGER, INTENT(IN) :: nspins, unit_nr, dimen
688 INTEGER, DIMENSION(nspins), INTENT(IN) :: virtual, homo
689 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: mo_coeff_o, mo_coeff_v
690 REAL(kind=dp), DIMENSION(dimen, nspins), &
691 INTENT(IN) :: eigenval
692 TYPE(qs_p_env_type), INTENT(IN), POINTER :: p_env
693 TYPE(cp_fm_type), DIMENSION(dimen), INTENT(IN) :: l_jb
694 TYPE(cp_fm_type), INTENT(INOUT) :: fm_g_mu_nu
695 TYPE(cp_fm_type), INTENT(IN) :: fm_back
696 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: p_ia
697
698 CHARACTER(LEN=*), PARAMETER :: routinen = 'solve_z_vector_eq_low'
699
700 INTEGER :: handle, i_global, iib, ispin, j_global, &
701 jjb, ncol_local, nrow_local
702 INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
703 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: precond
704
705 CALL timeset(routinen, handle)
706
707 ! Pople method
708 ! change sign to L_jb
709 DO ispin = 1, nspins
710 l_jb(ispin)%local_data(:, :) = -l_jb(ispin)%local_data(:, :)
711 END DO
712
713 ! create fm structure
714 ALLOCATE (precond(nspins))
715 DO ispin = 1, nspins
716 ! create preconditioner (for now only orbital energy differences)
717 CALL cp_fm_create(precond(ispin), p_ia(ispin)%matrix_struct, name="precond")
718 CALL cp_fm_set_all(precond(ispin), 1.0_dp)
719 CALL cp_fm_get_info(matrix=precond(ispin), &
720 nrow_local=nrow_local, &
721 ncol_local=ncol_local, &
722 row_indices=row_indices, &
723 col_indices=col_indices)
724 DO jjb = 1, ncol_local
725 j_global = col_indices(jjb)
726 DO iib = 1, nrow_local
727 i_global = row_indices(iib)
728 precond(ispin)%local_data(iib, jjb) = eigenval(j_global + homo(ispin), ispin) - eigenval(i_global, ispin)
729 END DO
730 END DO
731 END DO
732
733 DO ispin = 1, nspins
734 precond(ispin)%local_data(:, :) = 1.0_dp/precond(ispin)%local_data(:, :)
735 END DO
736
737 SELECT CASE (mp2_env%ri_grad%z_solver_method)
738 CASE (z_solver_pople)
739 CALL solve_z_vector_pople(qs_env, mp2_env, homo, virtual, dimen, unit_nr, nspins, &
740 mo_coeff_o, mo_coeff_v, eigenval, p_env, &
741 l_jb, fm_g_mu_nu, fm_back, p_ia, precond)
742 CASE (z_solver_cg, z_solver_richardson, z_solver_sd)
743 CALL solve_z_vector_cg(qs_env, mp2_env, homo, virtual, dimen, unit_nr, nspins, &
744 mo_coeff_o, mo_coeff_v, eigenval, p_env, &
745 l_jb, fm_g_mu_nu, fm_back, p_ia, precond)
746 CASE DEFAULT
747 cpabort("Unknown solver")
748 END SELECT
749
750 CALL cp_fm_release(precond)
751
752 CALL timestop(handle)
753
754 END SUBROUTINE solve_z_vector_eq_low
755
756! **************************************************************************************************
757!> \brief ...
758!> \param qs_env ...
759!> \param mp2_env ...
760!> \param homo ...
761!> \param virtual ...
762!> \param dimen ...
763!> \param unit_nr ...
764!> \param nspins ...
765!> \param mo_coeff_o ...
766!> \param mo_coeff_v ...
767!> \param Eigenval ...
768!> \param p_env ...
769!> \param L_jb ...
770!> \param fm_G_mu_nu ...
771!> \param fm_back ...
772!> \param P_ia ...
773!> \param precond ...
774! **************************************************************************************************
775 SUBROUTINE solve_z_vector_pople(qs_env, mp2_env, homo, virtual, dimen, unit_nr, nspins, &
776 mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
777 L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
778 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
779 TYPE(mp2_type), INTENT(IN) :: mp2_env
780 INTEGER, INTENT(IN) :: nspins, unit_nr, dimen
781 INTEGER, DIMENSION(nspins), INTENT(IN) :: virtual, homo
782 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: mo_coeff_o, mo_coeff_v
783 REAL(kind=dp), DIMENSION(dimen, nspins), &
784 INTENT(IN) :: eigenval
785 TYPE(qs_p_env_type), INTENT(IN), POINTER :: p_env
786 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: l_jb
787 TYPE(cp_fm_type), INTENT(INOUT) :: fm_g_mu_nu
788 TYPE(cp_fm_type), INTENT(IN) :: fm_back
789 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: p_ia, precond
790
791 CHARACTER(LEN=*), PARAMETER :: routinen = 'solve_z_vector_pople'
792
793 INTEGER :: cycle_counter, handle, iib, iiter, &
794 ispin, max_num_iter, transf_type_in
795 LOGICAL :: converged
796 REAL(kind=dp) :: conv, eps_conv, scale_cphf, t1, t2
797 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: proj_bi_xj, temp_vals, x_norm, xi_b
798 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: a_small, b_small, xi_axi
799 TYPE(cp_fm_struct_p_type), ALLOCATABLE, &
800 DIMENSION(:) :: fm_struct_tmp
801 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: b_i, residual
802 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: ax, xn
803
804 CALL timeset(routinen, handle)
805
806 eps_conv = mp2_env%ri_grad%cphf_eps_conv
807
808 IF (sqrt(accurate_dot_product_spin(l_jb, l_jb)) >= eps_conv) THEN
809
810 max_num_iter = mp2_env%ri_grad%cphf_max_num_iter
811 scale_cphf = mp2_env%ri_grad%scale_step_size
812
813 ! set the transformation type (equal for all methods all updates)
814 transf_type_in = 3
815
816 ! set convergence flag
817 converged = .false.
818
819 ALLOCATE (fm_struct_tmp(nspins), b_i(nspins), residual(nspins))
820 DO ispin = 1, nspins
821 fm_struct_tmp(ispin)%struct => p_ia(ispin)%matrix_struct
822
823 CALL cp_fm_create(b_i(ispin), fm_struct_tmp(ispin)%struct, name="b_i")
824 CALL cp_fm_set_all(b_i(ispin), 0.0_dp)
825 b_i(ispin)%local_data(:, :) = precond(ispin)%local_data(:, :)*l_jb(ispin)%local_data(:, :)
826
827 ! create the residual vector (r), we check convergence on the norm of
828 ! this vector r=(Ax-b)
829 CALL cp_fm_create(residual(ispin), fm_struct_tmp(ispin)%struct, name="residual")
830 CALL cp_fm_set_all(residual(ispin), 0.0_dp)
831 END DO
832
833 IF (unit_nr > 0) THEN
834 WRITE (unit_nr, *)
835 WRITE (unit_nr, '(T3,A)') "MP2_CPHF| Iterative solution of Z-Vector equations (Pople's method)"
836 WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', eps_conv
837 WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', max_num_iter
838 WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Scaling of initial guess: ', scale_cphf
839 WRITE (unit_nr, '(T4,A)') repeat("-", 40)
840 WRITE (unit_nr, '(T4,A,T15,A,T33,A)') 'Step', 'Time', 'Convergence'
841 WRITE (unit_nr, '(T4,A)') repeat("-", 40)
842 END IF
843
844 ALLOCATE (xn(nspins, max_num_iter))
845 ALLOCATE (ax(nspins, max_num_iter))
846 ALLOCATE (x_norm(max_num_iter))
847 ALLOCATE (xi_b(max_num_iter))
848 ALLOCATE (xi_axi(max_num_iter, 0:max_num_iter))
849 x_norm = 0.0_dp
850 xi_b = 0.0_dp
851 xi_axi = 0.0_dp
852
853 cycle_counter = 0
854 DO iiter = 1, max_num_iter
855 cycle_counter = cycle_counter + 1
856
857 t1 = m_walltime()
858
859 ! create and update x_i (orthogonalization with previous vectors)
860 DO ispin = 1, nspins
861 CALL cp_fm_create(xn(ispin, iiter), fm_struct_tmp(ispin)%struct, name="xi")
862 CALL cp_fm_set_all(xn(ispin, iiter), 0.0_dp)
863 END DO
864
865 ALLOCATE (proj_bi_xj(iiter - 1))
866 proj_bi_xj = 0.0_dp
867 ! first compute the projection of the actual b_i into all previous x_i
868 ! already scaled with the norm of each x_i
869 DO iib = 1, iiter - 1
870 proj_bi_xj(iib) = proj_bi_xj(iib) + accurate_dot_product_spin(b_i, xn(:, iib))/x_norm(iib)
871 END DO
872
873 ! update actual x_i
874 DO ispin = 1, nspins
875 xn(ispin, iiter)%local_data(:, :) = scale_cphf*b_i(ispin)%local_data(:, :)
876 DO iib = 1, iiter - 1
877 xn(ispin, iiter)%local_data(:, :) = xn(ispin, iiter)%local_data(:, :) - &
878 xn(ispin, iib)%local_data(:, :)*proj_bi_xj(iib)
879 END DO
880 END DO
881 DEALLOCATE (proj_bi_xj)
882
883 ! create Ax(iiter) that will store the matrix vector product for this cycle
884 DO ispin = 1, nspins
885 CALL cp_fm_create(ax(ispin, iiter), fm_struct_tmp(ispin)%struct, name="Ai")
886 CALL cp_fm_set_all(ax(ispin, iiter), 0.0_dp)
887 END DO
888
889 CALL cphf_like_update(qs_env, homo, virtual, dimen, nspins, &
890 mo_coeff_o, &
891 mo_coeff_v, eigenval, p_env, &
892 xn(:, iiter), fm_g_mu_nu, fm_back, transf_type_in, &
893 ax(:, iiter))
894
895 ! in order to reduce the number of parallel sums here we
896 ! cluster all necessary scalar products into a single vector
897 ! temp_vals contains:
898 ! 1:iiter -> <Ax_i|x_j>
899 ! iiter+1 -> <x_i|b>
900 ! iiter+2 -> <x_i|x_i>
901
902 ALLOCATE (temp_vals(iiter + 2))
903 temp_vals = 0.0_dp
904 ! <Ax_i|x_j>
905 DO iib = 1, iiter
906 temp_vals(iib) = temp_vals(iib) + accurate_dot_product_spin(ax(:, iiter), xn(:, iib))
907 END DO
908 ! <x_i|b>
909 temp_vals(iiter + 1) = temp_vals(iiter + 1) + accurate_dot_product_spin(xn(:, iiter), l_jb)
910 ! norm
911 temp_vals(iiter + 2) = temp_vals(iiter + 2) + accurate_dot_product_spin(xn(:, iiter), xn(:, iiter))
912 ! update <Ax_i|x_j>, <x_i|b> and norm <x_i|x_i>
913 xi_axi(iiter, 1:iiter) = temp_vals(1:iiter)
914 xi_axi(1:iiter, iiter) = temp_vals(1:iiter)
915 xi_b(iiter) = temp_vals(iiter + 1)
916 x_norm(iiter) = temp_vals(iiter + 2)
917 DEALLOCATE (temp_vals)
918
919 ! solve reduced system
920 IF (ALLOCATED(a_small)) DEALLOCATE (a_small)
921 IF (ALLOCATED(b_small)) DEALLOCATE (b_small)
922 ALLOCATE (a_small(iiter, iiter))
923 ALLOCATE (b_small(iiter, 1))
924 a_small(1:iiter, 1:iiter) = xi_axi(1:iiter, 1:iiter)
925 b_small(1:iiter, 1) = xi_b(1:iiter)
926
927 CALL solve_system(matrix=a_small, mysize=iiter, eigenvectors=b_small)
928
929 ! check for convergence
930 DO ispin = 1, nspins
931 CALL cp_fm_set_all(residual(ispin), 0.0_dp)
932 DO iib = 1, iiter
933 residual(ispin)%local_data(:, :) = &
934 residual(ispin)%local_data(:, :) + &
935 b_small(iib, 1)*ax(ispin, iib)%local_data(:, :)
936 END DO
937
938 residual(ispin)%local_data(:, :) = &
939 residual(ispin)%local_data(:, :) - &
940 l_jb(ispin)%local_data(:, :)
941 END DO
942
943 conv = sqrt(accurate_dot_product_spin(residual, residual))
944
945 t2 = m_walltime()
946
947 IF (unit_nr > 0) THEN
948 WRITE (unit_nr, '(T3,I5,T13,F6.1,11X,F14.8)') iiter, t2 - t1, conv
949 CALL m_flush(unit_nr)
950 END IF
951
952 IF (conv <= eps_conv) THEN
953 converged = .true.
954 EXIT
955 END IF
956
957 ! update b_i for the next round
958 DO ispin = 1, nspins
959 b_i(ispin)%local_data(:, :) = b_i(ispin)%local_data(:, :) &
960 + precond(ispin)%local_data(:, :) &
961 *ax(ispin, iiter)%local_data(:, :)
962 END DO
963
964 scale_cphf = 1.0_dp
965
966 END DO
967
968 IF (unit_nr > 0) THEN
969 WRITE (unit_nr, '(T4,A)') repeat("-", 40)
970 IF (converged) THEN
971 WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations converged in', cycle_counter, ' steps'
972 ELSE
973 WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations NOT converged in', cycle_counter, ' steps'
974 END IF
975 END IF
976
977 ! store solution into P_ia
978 DO iiter = 1, cycle_counter
979 DO ispin = 1, nspins
980 p_ia(ispin)%local_data(:, :) = p_ia(ispin)%local_data(:, :) + &
981 b_small(iiter, 1)*xn(ispin, iiter)%local_data(:, :)
982 END DO
983 END DO
984
985 ! Release arrays
986 DEALLOCATE (x_norm)
987 DEALLOCATE (xi_b)
988 DEALLOCATE (xi_axi)
989
990 CALL cp_fm_release(b_i)
991 CALL cp_fm_release(residual)
992 CALL cp_fm_release(ax)
993 CALL cp_fm_release(xn)
994
995 ELSE
996 IF (unit_nr > 0) THEN
997 WRITE (unit_nr, '(T4,A)') repeat("-", 40)
998 WRITE (unit_nr, '(T3,A)') 'Residual smaller than EPS_CONV. Skip solution of Z-vector equation.'
999 END IF
1000 END IF
1001
1002 CALL timestop(handle)
1003
1004 END SUBROUTINE solve_z_vector_pople
1005
1006! **************************************************************************************************
1007!> \brief ...
1008!> \param qs_env ...
1009!> \param mp2_env ...
1010!> \param homo ...
1011!> \param virtual ...
1012!> \param dimen ...
1013!> \param unit_nr ...
1014!> \param nspins ...
1015!> \param mo_coeff_o ...
1016!> \param mo_coeff_v ...
1017!> \param Eigenval ...
1018!> \param p_env ...
1019!> \param L_jb ...
1020!> \param fm_G_mu_nu ...
1021!> \param fm_back ...
1022!> \param P_ia ...
1023!> \param precond ...
1024! **************************************************************************************************
1025 SUBROUTINE solve_z_vector_cg(qs_env, mp2_env, homo, virtual, dimen, unit_nr, nspins, &
1026 mo_coeff_o, mo_coeff_v, Eigenval, p_env, &
1027 L_jb, fm_G_mu_nu, fm_back, P_ia, precond)
1028 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1029 TYPE(mp2_type), INTENT(IN) :: mp2_env
1030 INTEGER, INTENT(IN) :: nspins, unit_nr, dimen
1031 INTEGER, DIMENSION(nspins), INTENT(IN) :: virtual, homo
1032 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: mo_coeff_o, mo_coeff_v
1033 REAL(kind=dp), DIMENSION(dimen, nspins), &
1034 INTENT(IN) :: eigenval
1035 TYPE(qs_p_env_type), INTENT(IN), POINTER :: p_env
1036 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: l_jb
1037 TYPE(cp_fm_type), INTENT(INOUT) :: fm_g_mu_nu
1038 TYPE(cp_fm_type), INTENT(IN) :: fm_back
1039 TYPE(cp_fm_type), DIMENSION(nspins), INTENT(IN) :: p_ia, precond
1040
1041 CHARACTER(LEN=*), PARAMETER :: routinen = 'solve_z_vector_cg'
1042
1043 INTEGER :: cycles_passed, handle, iiter, ispin, max_num_iter, restart_counter, &
1044 restart_every, transf_type_in, z_solver_method
1045 LOGICAL :: converged, do_restart
1046 REAL(kind=dp) :: eps_conv, norm_residual, norm_residual_old, &
1047 residual_dot_diff_search_vec_new, residual_dot_diff_search_vec_old, &
1048 residual_dot_search_vec, residual_new_dot_diff_search_vec_old, scale_result, &
1049 scale_search, scale_step_size, search_vec_dot_a_search_vec, t1, t2
1050 TYPE(cp_fm_struct_p_type), ALLOCATABLE, &
1051 DIMENSION(:) :: fm_struct_tmp
1052 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: a_dot_search_vector, diff_search_vector, &
1053 residual, search_vector
1054
1055 CALL timeset(routinen, handle)
1056
1057 max_num_iter = mp2_env%ri_grad%cphf_max_num_iter
1058 eps_conv = mp2_env%ri_grad%cphf_eps_conv
1059 z_solver_method = mp2_env%ri_grad%z_solver_method
1060 restart_every = mp2_env%ri_grad%cphf_restart
1061 scale_step_size = mp2_env%ri_grad%scale_step_size
1062 transf_type_in = 3
1063
1064 IF (unit_nr > 0) THEN
1065 WRITE (unit_nr, *)
1066 SELECT CASE (z_solver_method)
1067 CASE (z_solver_cg)
1068 IF (mp2_env%ri_grad%polak_ribiere) THEN
1069 WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (CG with Polak-Ribiere step)'
1070 ELSE
1071 WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (CG with Fletcher-Reeves step)'
1072 END IF
1073 CASE (z_solver_richardson)
1074 WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (Richardson method)'
1075 CASE (z_solver_sd)
1076 WRITE (unit_nr, '(T3,A)') 'MP2_CPHF| Iterative solution of Z-Vector equations (Steepest Descent method)'
1077 CASE DEFAULT
1078 cpabort("Unknown solver")
1079 END SELECT
1080 WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Convergence threshold:', eps_conv
1081 WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Maximum number of iterations: ', max_num_iter
1082 WRITE (unit_nr, '(T3,A,T45,I8)') 'MP2_CPHF| Number of steps for restart: ', restart_every
1083 WRITE (unit_nr, '(T3, A)') 'MP2_CPHF| Restart after no decrease'
1084 WRITE (unit_nr, '(T3,A,T45,ES8.1)') 'MP2_CPHF| Scaling factor of each step: ', scale_step_size
1085 WRITE (unit_nr, '(T4,A)') repeat("-", 40)
1086 WRITE (unit_nr, '(T4,A,T13,A,T28,A,T43,A)') 'Step', 'Restart', 'Time', 'Convergence'
1087 WRITE (unit_nr, '(T4,A)') repeat("-", 40)
1088 END IF
1089
1090 ALLOCATE (fm_struct_tmp(nspins), residual(nspins), diff_search_vector(nspins), &
1091 search_vector(nspins), a_dot_search_vector(nspins))
1092 DO ispin = 1, nspins
1093 fm_struct_tmp(ispin)%struct => p_ia(ispin)%matrix_struct
1094
1095 CALL cp_fm_create(residual(ispin), fm_struct_tmp(ispin)%struct, name="residual")
1096 CALL cp_fm_set_all(residual(ispin), 0.0_dp)
1097
1098 CALL cp_fm_create(diff_search_vector(ispin), fm_struct_tmp(ispin)%struct, name="difference search vector")
1099 CALL cp_fm_set_all(diff_search_vector(ispin), 0.0_dp)
1100
1101 CALL cp_fm_create(search_vector(ispin), fm_struct_tmp(ispin)%struct, name="search vector")
1102 CALL cp_fm_set_all(search_vector(ispin), 0.0_dp)
1103
1104 CALL cp_fm_create(a_dot_search_vector(ispin), fm_struct_tmp(ispin)%struct, name="A times search vector")
1105 CALL cp_fm_set_all(a_dot_search_vector(ispin), 0.0_dp)
1106 END DO
1107
1108 converged = .false.
1109 cycles_passed = max_num_iter
1110 ! By that, we enforce the setup of the matrices
1111 do_restart = .true.
1112
1113 t1 = m_walltime()
1114
1115 DO iiter = 1, max_num_iter
1116
1117 ! During the first iteration, P_ia=0 such that the application of the 2nd order matrix is zero
1118 IF (do_restart) THEN
1119 ! We do not consider the first step to be a restart
1120 ! Do not recalculate residual if it is already enforced to save FLOPs
1121 IF (.NOT. mp2_env%ri_grad%recalc_residual .OR. (iiter == 1)) THEN
1122 IF (iiter > 1) THEN
1123 CALL cphf_like_update(qs_env, homo, virtual, dimen, nspins, &
1124 mo_coeff_o, &
1125 mo_coeff_v, eigenval, p_env, &
1126 p_ia, fm_g_mu_nu, fm_back, transf_type_in, &
1127 residual)
1128 ELSE
1129 do_restart = .false.
1130
1131 DO ispin = 1, nspins
1132 CALL cp_fm_set_all(residual(ispin), 0.0_dp)
1133 END DO
1134 END IF
1135
1136 DO ispin = 1, nspins
1137 residual(ispin)%local_data(:, :) = l_jb(ispin)%local_data(:, :) &
1138 - residual(ispin)%local_data(:, :)
1139 END DO
1140 END IF
1141
1142 DO ispin = 1, nspins
1143 diff_search_vector(ispin)%local_data(:, :) = &
1144 precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1145 search_vector(ispin)%local_data(:, :) = diff_search_vector(ispin)%local_data(:, :)
1146 END DO
1147
1148 restart_counter = 1
1149 END IF
1150
1151 norm_residual_old = sqrt(accurate_dot_product_spin(residual, residual))
1152
1153 residual_dot_diff_search_vec_old = accurate_dot_product_spin(residual, diff_search_vector)
1154
1155 CALL cphf_like_update(qs_env, homo, virtual, dimen, nspins, &
1156 mo_coeff_o, &
1157 mo_coeff_v, eigenval, p_env, &
1158 search_vector, fm_g_mu_nu, fm_back, transf_type_in, &
1159 a_dot_search_vector)
1160
1161 IF (z_solver_method /= z_solver_richardson) THEN
1162 search_vec_dot_a_search_vec = accurate_dot_product_spin(search_vector, a_dot_search_vector)
1163
1164 IF (z_solver_method == z_solver_cg) THEN
1165 scale_result = residual_dot_diff_search_vec_old/search_vec_dot_a_search_vec
1166 ELSE
1167 residual_dot_search_vec = accurate_dot_product_spin(residual, search_vector)
1168 scale_result = residual_dot_search_vec/search_vec_dot_a_search_vec
1169 END IF
1170
1171 scale_result = scale_result*scale_step_size
1172
1173 ELSE
1174
1175 scale_result = scale_step_size
1176
1177 END IF
1178
1179 DO ispin = 1, nspins
1180 p_ia(ispin)%local_data(:, :) = p_ia(ispin)%local_data(:, :) &
1181 + scale_result*search_vector(ispin)%local_data(:, :)
1182 END DO
1183
1184 IF (.NOT. mp2_env%ri_grad%recalc_residual) THEN
1185
1186 DO ispin = 1, nspins
1187 residual(ispin)%local_data(:, :) = residual(ispin)%local_data(:, :) &
1188 - scale_result*a_dot_search_vector(ispin)%local_data(:, :)
1189 END DO
1190 ELSE
1191 CALL cphf_like_update(qs_env, homo, virtual, dimen, nspins, &
1192 mo_coeff_o, &
1193 mo_coeff_v, eigenval, p_env, &
1194 p_ia, fm_g_mu_nu, fm_back, transf_type_in, &
1195 residual)
1196
1197 DO ispin = 1, nspins
1198 residual(ispin)%local_data(:, :) = l_jb(ispin)%local_data(:, :) - residual(ispin)%local_data(:, :)
1199 END DO
1200 END IF
1201
1202 norm_residual = sqrt(accurate_dot_product_spin(residual, residual))
1203
1204 t2 = m_walltime()
1205
1206 IF (unit_nr > 0) THEN
1207 WRITE (unit_nr, '(T3,I4,T16,L1,T26,F6.1,8X,F14.8)') iiter, do_restart, t2 - t1, norm_residual
1208 CALL m_flush(unit_nr)
1209 END IF
1210
1211 IF (norm_residual <= eps_conv) THEN
1212 converged = .true.
1213 cycles_passed = iiter
1214 EXIT
1215 END IF
1216
1217 t1 = m_walltime()
1218
1219 IF (z_solver_method == z_solver_richardson) THEN
1220 DO ispin = 1, nspins
1221 search_vector(ispin)%local_data(:, :) = &
1222 scale_step_size*precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1223 END DO
1224 ELSE IF (z_solver_method == z_solver_sd) THEN
1225 DO ispin = 1, nspins
1226 search_vector(ispin)%local_data(:, :) = &
1227 precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1228 END DO
1229 ELSE
1230 IF (mp2_env%ri_grad%polak_ribiere) &
1231 residual_new_dot_diff_search_vec_old = accurate_dot_product_spin(residual, diff_search_vector)
1232
1233 DO ispin = 1, nspins
1234 diff_search_vector(ispin)%local_data(:, :) = &
1235 precond(ispin)%local_data(:, :)*residual(ispin)%local_data(:, :)
1236 END DO
1237
1238 residual_dot_diff_search_vec_new = accurate_dot_product_spin(residual, diff_search_vector)
1239
1240 scale_search = residual_dot_diff_search_vec_new/residual_dot_diff_search_vec_old
1241 IF (mp2_env%ri_grad%polak_ribiere) scale_search = &
1242 scale_search - residual_new_dot_diff_search_vec_old/residual_dot_diff_search_vec_old
1243
1244 DO ispin = 1, nspins
1245 search_vector(ispin)%local_data(:, :) = scale_search*search_vector(ispin)%local_data(:, :) &
1246 + diff_search_vector(ispin)%local_data(:, :)
1247 END DO
1248
1249 ! Make new to old
1250 residual_dot_diff_search_vec_old = residual_dot_diff_search_vec_new
1251 END IF
1252
1253 ! Check whether the residual decrease or restart is enforced and ask for restart
1254 do_restart = (norm_residual >= norm_residual_old .OR. (mod(restart_counter, restart_every) == 0))
1255
1256 restart_counter = restart_counter + 1
1257 norm_residual_old = norm_residual
1258
1259 END DO
1260
1261 IF (unit_nr > 0) THEN
1262 WRITE (unit_nr, '(T4,A)') repeat("-", 40)
1263 IF (converged) THEN
1264 WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations converged in', cycles_passed, ' steps'
1265 ELSE
1266 WRITE (unit_nr, '(T3,A,I5,A)') 'Z-Vector equations NOT converged in', max_num_iter, ' steps'
1267 END IF
1268 END IF
1269
1270 DEALLOCATE (fm_struct_tmp)
1271 CALL cp_fm_release(residual)
1272 CALL cp_fm_release(diff_search_vector)
1273 CALL cp_fm_release(search_vector)
1274 CALL cp_fm_release(a_dot_search_vector)
1275
1276 CALL timestop(handle)
1277
1278 END SUBROUTINE solve_z_vector_cg
1279
1280! **************************************************************************************************
1281!> \brief ...
1282!> \param matrix1 ...
1283!> \param matrix2 ...
1284!> \return ...
1285! **************************************************************************************************
1286 FUNCTION accurate_dot_product_spin(matrix1, matrix2) RESULT(dotproduct)
1287 TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: matrix1, matrix2
1288 REAL(kind=dp) :: dotproduct
1289
1290 INTEGER :: ispin
1291
1292 dotproduct = 0.0_dp
1293 DO ispin = 1, SIZE(matrix1)
1294 dotproduct = dotproduct + accurate_dot_product(matrix1(ispin)%local_data, matrix2(ispin)%local_data)
1295 END DO
1296 CALL matrix1(1)%matrix_struct%para_env%sum(dotproduct)
1297
1298 END FUNCTION accurate_dot_product_spin
1299
1300! **************************************************************************************************
1301!> \brief ...
1302!> \param qs_env ...
1303! **************************************************************************************************
1304 SUBROUTINE update_mp2_forces(qs_env)
1305 TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
1306
1307 CHARACTER(LEN=*), PARAMETER :: routinen = 'update_mp2_forces'
1308
1309 INTEGER :: alpha, beta, handle, idir, iounit, &
1310 ispin, nimages, nocc, nspins
1311 INTEGER, DIMENSION(3) :: comp
1312 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1313 LOGICAL :: do_exx, do_hfx, use_virial
1314 REAL(kind=dp) :: e_dummy, e_hartree, e_xc, ehartree, exc
1315 REAL(kind=dp), DIMENSION(3) :: deb
1316 REAL(kind=dp), DIMENSION(3, 3) :: h_stress, pv_virial
1317 TYPE(admm_type), POINTER :: admm_env
1318 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1319 TYPE(cell_type), POINTER :: cell
1320 TYPE(cp_logger_type), POINTER :: logger
1321 TYPE(dbcsr_p_type), ALLOCATABLE, DIMENSION(:), &
1322 TARGET :: matrix_ks_aux
1323 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_p_mp2, &
1324 matrix_p_mp2_admm, matrix_s, rho1, &
1325 rho_ao, rho_ao_aux, scrm
1326 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao_kp, scrm_kp
1327 TYPE(dft_control_type), POINTER :: dft_control
1328 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
1329 TYPE(linres_control_type), POINTER :: linres_control
1330 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
1331 TYPE(mp_para_env_type), POINTER :: para_env
1332 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1333 POINTER :: sab_orb, sac_ae, sac_ppl, sap_ppnl
1334 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1335 TYPE(pw_c1d_gs_type) :: pot_g, rho_tot_g, temp_pw_g
1336 TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: dvg
1337 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rho_mp2_g
1338 TYPE(pw_c1d_gs_type), POINTER :: rho_core
1339 TYPE(pw_env_type), POINTER :: pw_env
1340 TYPE(pw_poisson_type), POINTER :: poisson_env
1341 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1342 TYPE(pw_r3d_rs_type) :: pot_r, vh_rspace, vhxc_rspace
1343 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_mp2_r, rho_mp2_r_aux, rho_r, &
1344 tau_mp2_r, vadmm_rspace, vtau_rspace, &
1345 vxc_rspace
1346 TYPE(qs_dispersion_type), POINTER :: dispersion_env
1347 TYPE(qs_energy_type), POINTER :: energy
1348 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1349 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1350 TYPE(qs_ks_env_type), POINTER :: ks_env
1351 TYPE(qs_p_env_type), POINTER :: p_env
1352 TYPE(qs_rho_type), POINTER :: rho, rho_aux
1353 TYPE(section_vals_type), POINTER :: hfx_section, hfx_sections, input, &
1354 xc_section
1355 TYPE(task_list_type), POINTER :: task_list_aux_fit
1356 TYPE(virial_type), POINTER :: virial
1357
1358 CALL timeset(routinen, handle)
1359
1360 NULLIFY (input, pw_env, matrix_s, rho, energy, force, virial, &
1361 matrix_p_mp2, matrix_p_mp2_admm, matrix_ks, rho_core)
1362 CALL get_qs_env(qs_env, &
1363 ks_env=ks_env, &
1364 dft_control=dft_control, &
1365 pw_env=pw_env, &
1366 input=input, &
1367 mos=mos, &
1368 para_env=para_env, &
1369 matrix_s=matrix_s, &
1370 matrix_ks=matrix_ks, &
1371 matrix_p_mp2=matrix_p_mp2, &
1372 matrix_p_mp2_admm=matrix_p_mp2_admm, &
1373 rho=rho, &
1374 cell=cell, &
1375 force=force, &
1376 virial=virial, &
1377 sab_orb=sab_orb, &
1378 energy=energy, &
1379 rho_core=rho_core, &
1380 x_data=x_data)
1381
1382 logger => cp_get_default_logger()
1383 iounit = cp_print_key_unit_nr(logger, input, "DFT%XC%WF_CORRELATION%PRINT", &
1384 extension=".mp2Log")
1385
1386 do_exx = .false.
1387 IF (qs_env%mp2_env%method == ri_rpa_method_gpw) THEN
1388 hfx_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC%WF_CORRELATION%RI_RPA%HF")
1389 CALL section_vals_get(hfx_section, explicit=do_exx)
1390 END IF
1391
1392 nimages = dft_control%nimages
1393 cpassert(nimages == 1)
1394 NULLIFY (cell_to_index)
1395
1396 p_env => qs_env%mp2_env%ri_grad%p_env
1397
1398 CALL qs_rho_get(rho, rho_ao=rho_ao, rho_ao_kp=rho_ao_kp, rho_r=rho_r, rho_g=rho_g)
1399 nspins = SIZE(rho_ao)
1400
1401 ! check if we have to calculate the virial
1402 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1403 IF (use_virial) virial%pv_calculate = .true.
1404
1405 CALL zero_qs_force(force)
1406 IF (use_virial) CALL zero_virial(virial, .false.)
1407
1408 DO ispin = 1, nspins
1409 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
1410 END DO
1411
1412 IF (nspins == 2) THEN
1413 CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, 1.0_dp, 1.0_dp)
1414 END IF
1415
1416 ! Kinetic energy matrix
1417 NULLIFY (scrm)
1418 IF (debug_forces) THEN
1419 deb(1:3) = force(1)%kinetic(1:3, 1)
1420 IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1421 END IF
1422 CALL build_kinetic_matrix(ks_env, matrix_t=scrm, &
1423 matrix_name="KINETIC ENERGY MATRIX", &
1424 basis_type="ORB", &
1425 sab_nl=sab_orb, calculate_forces=.true., &
1426 matrix_p=rho_ao(1)%matrix)
1427 IF (debug_forces) THEN
1428 deb(1:3) = force(1)%kinetic(1:3, 1) - deb(1:3)
1429 CALL para_env%sum(deb)
1430 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dT ", deb
1431 IF (use_virial) THEN
1432 e_dummy = third_tr(virial%pv_virial) - e_dummy
1433 CALL para_env%sum(e_dummy)
1434 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: T ", e_dummy
1435 END IF
1436 END IF
1437 CALL dbcsr_deallocate_matrix_set(scrm)
1438
1439 IF (nspins == 2) THEN
1440 CALL dbcsr_add(rho_ao(1)%matrix, rho_ao(2)%matrix, 1.0_dp, -1.0_dp)
1441 END IF
1442
1443 ! Add pseudo potential terms
1444 scrm_kp(1:nspins, 1:1) => matrix_ks(1:nspins)
1445 CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, particle_set=particle_set, &
1446 atomic_kind_set=atomic_kind_set, sac_ae=sac_ae, sac_ppl=sac_ppl, sap_ppnl=sap_ppnl)
1447 IF (ASSOCIATED(sac_ae)) THEN
1448 IF (debug_forces) THEN
1449 deb(1:3) = force(1)%all_potential(1:3, 1)
1450 IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1451 END IF
1452 CALL build_core_ae(scrm_kp, rho_ao_kp, force, &
1453 virial, .true., use_virial, 1, &
1454 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ae, &
1455 nimages, cell_to_index)
1456 IF (debug_forces) THEN
1457 deb(1:3) = force(1)%all_potential(1:3, 1) - deb(1:3)
1458 CALL para_env%sum(deb)
1459 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHae ", deb
1460 IF (use_virial) THEN
1461 e_dummy = third_tr(virial%pv_virial) - e_dummy
1462 CALL para_env%sum(e_dummy)
1463 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hae ", e_dummy
1464 END IF
1465 END IF
1466 END IF
1467 IF (ASSOCIATED(sac_ppl)) THEN
1468 IF (debug_forces) THEN
1469 deb(1:3) = force(1)%gth_ppl(1:3, 1)
1470 IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1471 END IF
1472 CALL build_core_ppl(scrm_kp, rho_ao_kp, force, &
1473 virial, .true., use_virial, 1, &
1474 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sac_ppl, &
1475 nimages=1, basis_type="ORB")
1476 IF (debug_forces) THEN
1477 deb(1:3) = force(1)%gth_ppl(1:3, 1) - deb(1:3)
1478 CALL para_env%sum(deb)
1479 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppl ", deb
1480 IF (use_virial) THEN
1481 e_dummy = third_tr(virial%pv_virial) - e_dummy
1482 CALL para_env%sum(e_dummy)
1483 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hppl ", e_dummy
1484 END IF
1485 END IF
1486 END IF
1487 IF (ASSOCIATED(sap_ppnl)) THEN
1488 IF (debug_forces) THEN
1489 deb(1:3) = force(1)%gth_ppnl(1:3, 1)
1490 IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1491 END IF
1492 CALL build_core_ppnl(scrm_kp, rho_ao_kp, force, &
1493 virial, .true., use_virial, 1, &
1494 qs_kind_set, atomic_kind_set, particle_set, sab_orb, sap_ppnl, dft_control%qs_control%eps_ppnl, &
1495 nimages=1, basis_type="ORB")
1496 IF (debug_forces) THEN
1497 deb(1:3) = force(1)%gth_ppnl(1:3, 1) - deb(1:3)
1498 CALL para_env%sum(deb)
1499 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dHppnl ", deb
1500 IF (use_virial) THEN
1501 e_dummy = third_tr(virial%pv_virial) - e_dummy
1502 CALL para_env%sum(e_dummy)
1503 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hppnl ", e_dummy
1504 END IF
1505 END IF
1506 END IF
1507
1508 ! Get the different components of the KS potential
1509 NULLIFY (vxc_rspace, vtau_rspace, vadmm_rspace)
1510 IF (use_virial) THEN
1511 h_stress = 0.0_dp
1512 CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc, h_stress)
1513 ! Update virial
1514 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1515 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1516 IF (.NOT. do_exx) THEN
1517 virial%pv_exc = virial%pv_exc - virial%pv_xc
1518 virial%pv_virial = virial%pv_virial - virial%pv_xc
1519 ELSE
1520 virial%pv_xc = 0.0_dp
1521 END IF
1522 IF (debug_forces) THEN
1523 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: xc ", third_tr(h_stress)
1524 CALL para_env%sum(virial%pv_xc(1, 1))
1525 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Corr xc ", third_tr(virial%pv_xc)
1526 END IF
1527 ELSE
1528 CALL ks_ref_potential(qs_env, vh_rspace, vxc_rspace, vtau_rspace, vadmm_rspace, ehartree, exc)
1529 END IF
1530
1531 ! Vhxc
1532 CALL get_qs_env(qs_env, pw_env=pw_env)
1533 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1534 poisson_env=poisson_env)
1535 CALL auxbas_pw_pool%create_pw(vhxc_rspace)
1536 IF (use_virial) h_stress = virial%pv_virial
1537 IF (debug_forces) THEN
1538 deb(1:3) = force(1)%rho_elec(1:3, 1)
1539 IF (use_virial) e_dummy = third_tr(h_stress)
1540 END IF
1541 IF (do_exx) THEN
1542 DO ispin = 1, nspins
1543 CALL pw_transfer(vh_rspace, vhxc_rspace)
1544 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
1545 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1546 hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
1547 qs_env=qs_env, calculate_forces=.true.)
1548 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, 1.0_dp)
1549 CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1550 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1551 hmat=matrix_ks(ispin), pmat=matrix_p_mp2(ispin), &
1552 qs_env=qs_env, calculate_forces=.true.)
1553 IF (ASSOCIATED(vtau_rspace)) THEN
1554 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
1555 hmat=matrix_ks(ispin), pmat=matrix_p_mp2(ispin), &
1556 qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
1557 END IF
1558 END DO
1559 ELSE
1560 DO ispin = 1, nspins
1561 CALL pw_transfer(vh_rspace, vhxc_rspace)
1562 CALL pw_axpy(vxc_rspace(ispin), vhxc_rspace)
1563 CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
1564 hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
1565 qs_env=qs_env, calculate_forces=.true.)
1566 IF (ASSOCIATED(vtau_rspace)) THEN
1567 CALL integrate_v_rspace(v_rspace=vtau_rspace(ispin), &
1568 hmat=matrix_ks(ispin), pmat=rho_ao(ispin), &
1569 qs_env=qs_env, calculate_forces=.true., compute_tau=.true.)
1570 END IF
1571 END DO
1572 END IF
1573 IF (debug_forces) THEN
1574 deb(1:3) = force(1)%rho_elec(1:3, 1) - deb(1:3)
1575 CALL para_env%sum(deb)
1576 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*dVhxc ", deb
1577 IF (use_virial) THEN
1578 e_dummy = third_tr(virial%pv_virial) - e_dummy
1579 CALL para_env%sum(e_dummy)
1580 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Vhxc ", e_dummy
1581 END IF
1582 END IF
1583 IF (use_virial) THEN
1584 h_stress = virial%pv_virial - h_stress
1585 virial%pv_ehartree = virial%pv_ehartree + h_stress
1586
1587 CALL qs_rho_get(p_env%rho1, rho_r=rho_mp2_r, tau_r=tau_mp2_r)
1588 e_xc = 0.0_dp
1589 DO ispin = 1, nspins
1590 ! The potentials have been scaled in ks_ref_potential, but for pw_integral_ab, we need the unscaled potentials
1591 CALL pw_scale(vxc_rspace(ispin), 1.0_dp/vxc_rspace(ispin)%pw_grid%dvol)
1592 e_xc = e_xc + pw_integral_ab(rho_mp2_r(ispin), vxc_rspace(ispin))
1593 IF (ASSOCIATED(vtau_rspace)) CALL pw_scale(vtau_rspace(ispin), 1.0_dp/vtau_rspace(ispin)%pw_grid%dvol)
1594 IF (ASSOCIATED(vtau_rspace)) e_xc = e_xc + pw_integral_ab(tau_mp2_r(ispin), vtau_rspace(ispin))
1595 END DO
1596 IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG VIRIAL:: vxc*d1 ", e_xc
1597 DO alpha = 1, 3
1598 virial%pv_exc(alpha, alpha) = virial%pv_exc(alpha, alpha) - e_xc/real(para_env%num_pe, dp)
1599 virial%pv_virial(alpha, alpha) = virial%pv_virial(alpha, alpha) - e_xc/real(para_env%num_pe, dp)
1600 END DO
1601 END IF
1602 DO ispin = 1, nspins
1603 CALL auxbas_pw_pool%give_back_pw(vxc_rspace(ispin))
1604 IF (ASSOCIATED(vtau_rspace)) THEN
1605 CALL auxbas_pw_pool%give_back_pw(vtau_rspace(ispin))
1606 END IF
1607 END DO
1608 DEALLOCATE (vxc_rspace)
1609 CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
1610 IF (ASSOCIATED(vtau_rspace)) DEALLOCATE (vtau_rspace)
1611
1612 DO ispin = 1, nspins
1613 CALL dbcsr_add(rho_ao(ispin)%matrix, matrix_p_mp2(ispin)%matrix, 1.0_dp, -1.0_dp)
1614 END DO
1615
1616 ! pw stuff
1617 NULLIFY (poisson_env, auxbas_pw_pool)
1618 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
1619 poisson_env=poisson_env)
1620
1621 ! get some of the grids ready
1622 CALL auxbas_pw_pool%create_pw(pot_r)
1623 CALL auxbas_pw_pool%create_pw(pot_g)
1624 CALL auxbas_pw_pool%create_pw(rho_tot_g)
1625
1626 CALL pw_zero(rho_tot_g)
1627
1628 CALL qs_rho_get(p_env%rho1, rho_r=rho_mp2_r, rho_g=rho_mp2_g)
1629 DO ispin = 1, nspins
1630 CALL pw_axpy(rho_mp2_g(ispin), rho_tot_g)
1631 END DO
1632
1633 IF (use_virial) THEN
1634 ALLOCATE (dvg(3))
1635 DO idir = 1, 3
1636 CALL auxbas_pw_pool%create_pw(dvg(idir))
1637 END DO
1638 CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g, dvhartree=dvg)
1639 ELSE
1640 CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g)
1641 END IF
1642
1643 CALL pw_transfer(pot_g, pot_r)
1644 CALL pw_scale(pot_r, pot_r%pw_grid%dvol)
1645 CALL pw_axpy(pot_r, vh_rspace)
1646
1647 ! calculate core forces
1648 CALL integrate_v_core_rspace(vh_rspace, qs_env)
1649 IF (debug_forces) THEN
1650 deb(:) = force(1)%rho_core(:, 1)
1651 CALL para_env%sum(deb)
1652 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: core ", deb
1653 IF (use_virial) THEN
1654 e_dummy = third_tr(virial%pv_virial) - e_dummy
1655 CALL para_env%sum(e_dummy)
1656 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: core ", e_dummy
1657 END IF
1658 END IF
1659 CALL auxbas_pw_pool%give_back_pw(vh_rspace)
1660
1661 IF (use_virial) THEN
1662 ! update virial if necessary with the volume term
1663 ! first create pw auxiliary stuff
1664 CALL auxbas_pw_pool%create_pw(temp_pw_g)
1665
1666 ! make a copy of the MP2 density in G space
1667 CALL pw_copy(rho_tot_g, temp_pw_g)
1668
1669 ! calculate total SCF density and potential
1670 CALL pw_copy(rho_g(1), rho_tot_g)
1671 IF (nspins == 2) CALL pw_axpy(rho_g(2), rho_tot_g)
1672 CALL pw_axpy(rho_core, rho_tot_g)
1673 CALL pw_poisson_solve(poisson_env, rho_tot_g, vhartree=pot_g)
1674
1675 ! finally update virial with the volume contribution
1676 e_hartree = pw_integral_ab(temp_pw_g, pot_g)
1677 IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: vh1*d0 ", e_hartree
1678
1679 h_stress = 0.0_dp
1680 DO alpha = 1, 3
1681 comp = 0
1682 comp(alpha) = 1
1683 CALL pw_copy(pot_g, rho_tot_g)
1684 CALL pw_derive(rho_tot_g, comp)
1685 h_stress(alpha, alpha) = -e_hartree
1686 DO beta = alpha, 3
1687 h_stress(alpha, beta) = h_stress(alpha, beta) &
1688 - 2.0_dp*pw_integral_ab(rho_tot_g, dvg(beta))/fourpi
1689 h_stress(beta, alpha) = h_stress(alpha, beta)
1690 END DO
1691 END DO
1692 IF (debug_forces .AND. iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Hartree ", third_tr(h_stress)
1693
1694 ! free stuff
1695 CALL auxbas_pw_pool%give_back_pw(temp_pw_g)
1696 DO idir = 1, 3
1697 CALL auxbas_pw_pool%give_back_pw(dvg(idir))
1698 END DO
1699 DEALLOCATE (dvg)
1700
1701 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
1702 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
1703
1704 END IF
1705
1706 DO ispin = 1, nspins
1707 CALL dbcsr_set(p_env%kpp1(ispin)%matrix, 0.0_dp)
1708 IF (dft_control%do_admm) CALL dbcsr_set(p_env%kpp1_admm(ispin)%matrix, 0.0_dp)
1709 END DO
1710
1711 CALL get_qs_env(qs_env=qs_env, linres_control=linres_control)
1712
1713 IF (dft_control%do_admm) THEN
1714 CALL get_qs_env(qs_env, admm_env=admm_env)
1715 xc_section => admm_env%xc_section_primary
1716 ELSE
1717 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1718 END IF
1719
1720 IF (use_virial) THEN
1721 h_stress = 0.0_dp
1722 pv_virial = virial%pv_virial
1723 END IF
1724 IF (debug_forces) THEN
1725 deb = force(1)%rho_elec(1:3, 1)
1726 IF (use_virial) e_dummy = third_tr(pv_virial)
1727 END IF
1728 CALL apply_2nd_order_kernel(qs_env, p_env, .false., .true., use_virial, h_stress)
1729 IF (use_virial) THEN
1730 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_virial)
1731 IF (debug_forces) THEN
1732 e_dummy = third_tr(virial%pv_virial - pv_virial)
1733 CALL para_env%sum(e_dummy)
1734 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Kh ", e_dummy
1735 END IF
1736 virial%pv_exc = virial%pv_exc + h_stress
1737 virial%pv_virial = virial%pv_virial + h_stress
1738 IF (debug_forces) THEN
1739 e_dummy = third_tr(h_stress)
1740 CALL para_env%sum(e_dummy)
1741 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Kxc ", e_dummy
1742 END IF
1743 END IF
1744 IF (debug_forces) THEN
1745 deb(:) = force(1)%rho_elec(:, 1) - deb
1746 CALL para_env%sum(deb)
1747 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P0*Khxc ", deb
1748 IF (use_virial) THEN
1749 e_dummy = third_tr(virial%pv_virial) - e_dummy
1750 CALL para_env%sum(e_dummy)
1751 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: Khxc ", e_dummy
1752 END IF
1753 END IF
1754
1755 ! hfx section
1756 NULLIFY (hfx_sections)
1757 hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
1758 CALL section_vals_get(hfx_sections, explicit=do_hfx)
1759 IF (do_hfx) THEN
1760 IF (do_exx) THEN
1761 IF (dft_control%do_admm) THEN
1762 CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
1763 CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux, rho_ao_kp=rho_ao_kp)
1764 rho1 => p_env%p1_admm
1765 ELSE
1766 rho1 => p_env%p1
1767 END IF
1768 ELSE
1769 IF (dft_control%do_admm) THEN
1770 CALL get_admm_env(qs_env%admm_env, rho_aux_fit=rho_aux)
1771 CALL qs_rho_get(rho_aux, rho_ao=rho_ao_aux, rho_ao_kp=rho_ao_kp)
1772 DO ispin = 1, nspins
1773 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
1774 END DO
1775 rho1 => p_env%p1_admm
1776 ELSE
1777 DO ispin = 1, nspins
1778 CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
1779 END DO
1780 rho1 => p_env%p1
1781 END IF
1782 END IF
1783
1784 IF (x_data(1, 1)%do_hfx_ri) THEN
1785
1786 CALL hfx_ri_update_forces(qs_env, x_data(1, 1)%ri_data, nspins, &
1787 x_data(1, 1)%general_parameter%fraction, &
1788 rho_ao=rho_ao_kp, rho_ao_resp=rho1, &
1789 use_virial=use_virial, resp_only=do_exx)
1790
1791 ELSE
1792 CALL derivatives_four_center(qs_env, rho_ao_kp, rho1, hfx_sections, para_env, &
1793 1, use_virial, resp_only=do_exx)
1794 END IF
1795
1796 IF (use_virial) THEN
1797 virial%pv_exx = virial%pv_exx - virial%pv_fock_4c
1798 virial%pv_virial = virial%pv_virial - virial%pv_fock_4c
1799 END IF
1800 IF (debug_forces) THEN
1801 deb(1:3) = force(1)%fock_4c(1:3, 1) - deb(1:3)
1802 CALL para_env%sum(deb)
1803 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*hfx ", deb
1804 IF (use_virial) THEN
1805 e_dummy = third_tr(virial%pv_fock_4c)
1806 CALL para_env%sum(e_dummy)
1807 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: hfx ", e_dummy
1808 END IF
1809 END IF
1810
1811 IF (.NOT. do_exx) THEN
1812 IF (dft_control%do_admm) THEN
1813 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1814 DO ispin = 1, nspins
1815 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
1816 END DO
1817 ELSE
1818 DO ispin = 1, nspins
1819 CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, -1.0_dp)
1820 END DO
1821 END IF
1822 END IF
1823
1824 IF (dft_control%do_admm) THEN
1825 IF (debug_forces) THEN
1826 deb = force(1)%overlap_admm(1:3, 1)
1827 IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1828 END IF
1829 ! The 2nd order kernel contains a factor of two in apply_xc_admm_ao which we don't need for the projection derivatives
1830 IF (nspins == 1) CALL dbcsr_scale(p_env%kpp1_admm(1)%matrix, 0.5_dp)
1831 CALL admm_projection_derivative(qs_env, p_env%kpp1_admm, rho_ao)
1832 IF (debug_forces) THEN
1833 deb(:) = force(1)%overlap_admm(:, 1) - deb
1834 CALL para_env%sum(deb)
1835 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*KADMM*dS'", deb
1836 IF (use_virial) THEN
1837 e_dummy = third_tr(virial%pv_virial) - e_dummy
1838 CALL para_env%sum(e_dummy)
1839 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: KADMM*S' ", e_dummy
1840 END IF
1841 END IF
1842
1843 ALLOCATE (matrix_ks_aux(nspins))
1844 DO ispin = 1, nspins
1845 NULLIFY (matrix_ks_aux(ispin)%matrix)
1846 ALLOCATE (matrix_ks_aux(ispin)%matrix)
1847 CALL dbcsr_copy(matrix_ks_aux(ispin)%matrix, p_env%kpp1_admm(ispin)%matrix)
1848 CALL dbcsr_set(matrix_ks_aux(ispin)%matrix, 0.0_dp)
1849 END DO
1850
1851 ! Calculate kernel
1852 CALL tddft_hfx_matrix(matrix_ks_aux, rho_ao_aux, qs_env, .false.)
1853
1854 IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
1855 CALL get_qs_env(qs_env, ks_env=ks_env)
1856 CALL get_admm_env(qs_env%admm_env, task_list_aux_fit=task_list_aux_fit)
1857
1858 DO ispin = 1, nspins
1859 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, 1.0_dp)
1860 END DO
1861
1862 IF (use_virial) THEN
1863 CALL qs_rho_get(p_env%rho1_admm, rho_r=rho_mp2_r_aux)
1864 e_xc = 0.0_dp
1865 DO ispin = 1, nspins
1866 e_xc = e_xc + pw_integral_ab(rho_mp2_r_aux(ispin), vadmm_rspace(ispin))
1867 END DO
1868
1869 e_xc = -e_xc/vadmm_rspace(1)%pw_grid%dvol/real(para_env%num_pe, dp)
1870
1871 ! Update the virial
1872 DO alpha = 1, 3
1873 virial%pv_exc(alpha, alpha) = virial%pv_exc(alpha, alpha) + e_xc
1874 virial%pv_virial(alpha, alpha) = virial%pv_virial(alpha, alpha) + e_xc
1875 END DO
1876 IF (debug_forces) THEN
1877 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: P1*VADMM ", e_xc
1878 END IF
1879 END IF
1880
1881 IF (use_virial) h_stress = virial%pv_virial
1882 IF (debug_forces) THEN
1883 deb = force(1)%rho_elec(1:3, 1)
1884 IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1885 END IF
1886 DO ispin = 1, nspins
1887 IF (do_exx) THEN
1888 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=matrix_ks_aux(ispin), qs_env=qs_env, &
1889 calculate_forces=.true., basis_type="AUX_FIT", &
1890 task_list_external=task_list_aux_fit, pmat=matrix_p_mp2_admm(ispin))
1891 ELSE
1892 CALL integrate_v_rspace(v_rspace=vadmm_rspace(ispin), hmat=matrix_ks_aux(ispin), qs_env=qs_env, &
1893 calculate_forces=.true., basis_type="AUX_FIT", &
1894 task_list_external=task_list_aux_fit, pmat=rho_ao_aux(ispin))
1895 END IF
1896 CALL auxbas_pw_pool%give_back_pw(vadmm_rspace(ispin))
1897 END DO
1898 IF (use_virial) virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - h_stress)
1899 DEALLOCATE (vadmm_rspace)
1900 IF (debug_forces) THEN
1901 deb(:) = force(1)%rho_elec(:, 1) - deb
1902 CALL para_env%sum(deb)
1903 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*VADMM' ", deb
1904 IF (use_virial) THEN
1905 e_dummy = third_tr(virial%pv_virial) - e_dummy
1906 CALL para_env%sum(e_dummy)
1907 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: VADMM' ", e_dummy
1908 END IF
1909 END IF
1910
1911 DO ispin = 1, nspins
1912 CALL dbcsr_add(rho_ao_aux(ispin)%matrix, p_env%p1_admm(ispin)%matrix, 1.0_dp, -1.0_dp)
1913 END DO
1914
1915 END IF
1916
1917 DO ispin = 1, nspins
1918 CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
1919 END DO
1920
1921 IF (debug_forces) THEN
1922 deb = force(1)%overlap_admm(1:3, 1)
1923 IF (use_virial) e_dummy = third_tr(virial%pv_virial)
1924 END IF
1925 ! Add the second half of the projector deriatives contracting the first order density matrix with the fockian in the auxiliary basis
1926 IF (do_exx) THEN
1927 CALL admm_projection_derivative(qs_env, matrix_ks_aux, matrix_p_mp2)
1928 ELSE
1929 CALL admm_projection_derivative(qs_env, matrix_ks_aux, rho_ao)
1930 END IF
1931 IF (debug_forces) THEN
1932 deb(:) = force(1)%overlap_admm(:, 1) - deb
1933 CALL para_env%sum(deb)
1934 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*VADMM*dS'", deb
1935 IF (use_virial) THEN
1936 e_dummy = third_tr(virial%pv_virial) - e_dummy
1937 CALL para_env%sum(e_dummy)
1938 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: VADMM*S' ", e_dummy
1939 END IF
1940 END IF
1941
1942 DO ispin = 1, nspins
1943 CALL dbcsr_add(rho_ao(ispin)%matrix, p_env%p1(ispin)%matrix, 1.0_dp, -1.0_dp)
1944 END DO
1945
1946 DO ispin = 1, nspins
1947 CALL dbcsr_release(matrix_ks_aux(ispin)%matrix)
1948 DEALLOCATE (matrix_ks_aux(ispin)%matrix)
1949 END DO
1950 DEALLOCATE (matrix_ks_aux)
1951 END IF
1952 END IF
1953
1954 CALL dbcsr_scale(p_env%w1(1)%matrix, -1.0_dp)
1955
1956 ! Finish matrix_w_mp2 with occ-occ block
1957 DO ispin = 1, nspins
1958 CALL get_mo_set(mo_set=mos(ispin), homo=nocc, nmo=alpha)
1959 CALL calculate_whz_matrix(mos(ispin)%mo_coeff, p_env%kpp1(ispin)%matrix, &
1960 p_env%w1(1)%matrix, 1.0_dp, nocc)
1961 END DO
1962
1963 IF (debug_forces .AND. use_virial) e_dummy = third_tr(virial%pv_virial)
1964
1965 NULLIFY (scrm)
1966 CALL build_overlap_matrix(ks_env, matrix_s=scrm, &
1967 matrix_name="OVERLAP MATRIX", &
1968 basis_type_a="ORB", basis_type_b="ORB", &
1969 sab_nl=sab_orb, calculate_forces=.true., &
1970 matrix_p=p_env%w1(1)%matrix)
1971 CALL dbcsr_deallocate_matrix_set(scrm)
1972
1973 IF (debug_forces) THEN
1974 deb = force(1)%overlap(1:3, 1)
1975 CALL para_env%sum(deb)
1976 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: W*dS ", deb
1977 IF (use_virial) THEN
1978 e_dummy = third_tr(virial%pv_virial) - e_dummy
1979 CALL para_env%sum(e_dummy)
1980 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: S ", e_dummy
1981 END IF
1982 END IF
1983
1984 CALL auxbas_pw_pool%give_back_pw(pot_r)
1985 CALL auxbas_pw_pool%give_back_pw(pot_g)
1986 CALL auxbas_pw_pool%give_back_pw(rho_tot_g)
1987
1988 ! Release linres stuff
1989 CALL p_env_release(p_env)
1990 DEALLOCATE (p_env)
1991 NULLIFY (qs_env%mp2_env%ri_grad%p_env)
1992
1993 CALL sum_qs_force(force, qs_env%mp2_env%ri_grad%mp2_force)
1994 CALL deallocate_qs_force(qs_env%mp2_env%ri_grad%mp2_force)
1995
1996 IF (use_virial) THEN
1997 virial%pv_mp2 = qs_env%mp2_env%ri_grad%mp2_virial
1998 virial%pv_virial = virial%pv_virial + qs_env%mp2_env%ri_grad%mp2_virial
1999 IF (debug_forces) THEN
2000 e_dummy = third_tr(virial%pv_mp2)
2001 CALL para_env%sum(e_dummy)
2002 IF (iounit > 0) WRITE (iounit, "(T3,A,T33,F16.8)") "DEBUG VIRIAL:: MP2nonsep ", e_dummy
2003 END IF
2004 END IF
2005 ! Rewind the change from the beginning
2006 IF (use_virial) virial%pv_calculate = .false.
2007
2008 ! Add the dispersion forces and virials
2009 CALL get_qs_env(qs_env, dispersion_env=dispersion_env)
2010 CALL calculate_dispersion_pairpot(qs_env, dispersion_env, e_dummy, .true.)
2011
2012 CALL cp_print_key_finished_output(iounit, logger, input, &
2013 "DFT%XC%WF_CORRELATION%PRINT")
2014
2015 CALL timestop(handle)
2016
2017 END SUBROUTINE update_mp2_forces
2018
2019! **************************************************************************************************
2020!> \brief Calculates the third of the trace of a 3x3 matrix, for debugging purposes
2021!> \param matrix ...
2022!> \return ...
2023! **************************************************************************************************
2024 PURE FUNCTION third_tr(matrix)
2025 REAL(kind=dp), DIMENSION(3, 3), INTENT(IN) :: matrix
2026 REAL(kind=dp) :: third_tr
2027
2028 third_tr = (matrix(1, 1) + matrix(2, 2) + matrix(3, 3))/3.0_dp
2029
2030 END FUNCTION third_tr
2031
2032END MODULE mp2_cphf
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:87
kahan_sum::accurate_dot_product
Definition kahan_sum.F:52
parallel_gemm_api::parallel_gemm
Definition parallel_gemm_api.F:38
pw_methods::pw_axpy
Definition pw_methods.F:164
pw_methods::pw_copy
Definition pw_methods.F:145
pw_methods::pw_integral_ab
Definition pw_methods.F:221
pw_methods::pw_scale
Definition pw_methods.F:93
pw_methods::pw_transfer
Definition pw_methods.F:303
pw_methods::pw_zero
Definition pw_methods.F:82
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
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, atcore)
...
Definition core_ae.F:86
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, atcore)
...
Definition core_ppl.F:97
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, atcore)
...
Definition core_ppnl.F:90
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1165
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:368
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1181
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1119
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:251
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:214
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:695
cp_fm_basic_linalg
basic linear algebra operations for full matrices
Definition cp_fm_basic_linalg.F:14
cp_fm_basic_linalg::cp_fm_uplo_to_full
subroutine, public cp_fm_uplo_to_full(matrix, work, uplo)
given a triangular matrix according to uplo, computes the corresponding full matrix
Definition cp_fm_basic_linalg.F:1747
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:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:351
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:1009
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:1472
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:528
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:164
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:853
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:1063
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, external_para_env)
Add the hfx contributions to the Hamiltonian.
Definition hfx_admm_utils.F:2478
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:3039
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:15
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:2522
hfx_types::alloc_containers
subroutine, public alloc_containers(data, bin_size)
...
Definition hfx_types.F:2905
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:1074
input_constants::z_solver_cg
integer, parameter, public z_solver_cg
Definition input_constants.F:1074
input_constants::ri_rpa_method_gpw
integer, parameter, public ri_rpa_method_gpw
Definition input_constants.F:1059
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition input_constants.F:869
input_constants::z_solver_sd
integer, parameter, public z_solver_sd
Definition input_constants.F:1074
input_constants::z_solver_richardson
integer, parameter, public z_solver_richardson
Definition input_constants.F:1074
input_constants::ot_precond_full_all
integer, parameter, public ot_precond_full_all
Definition input_constants.F:514
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:731
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:704
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:15
linear_systems::solve_system
subroutine, public solve_system(matrix, mysize, eigenvectors)
...
Definition linear_systems.F:37
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:130
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:147
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:1305
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:10256
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:422
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:372
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_pp, 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, harris_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, eeq, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:514
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, harris_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, eeq, rhs)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1055
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_pp, 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:572
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:479
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:130
qs_p_env_methods::p_env_update_rho
subroutine, public p_env_update_rho(p_env, qs_env)
...
Definition qs_p_env_methods.F:395
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:314
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
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_blacs_env::cp_blacs_env_type
represent a blacs multidimensional parallel environment (for the mpi corrispective see cp_paratypes/m...
Definition cp_blacs_env.F:53
cp_control_types::dft_control_type
Definition cp_control_types.F:570
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:170
cp_fm_struct::cp_fm_struct_p_type
Definition cp_fm_struct.F:100
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:113
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
hfx_types::hfx_container_type
Definition hfx_types.F:253
hfx_types::hfx_type
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:509
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
mp2_types::mp2_type
Definition mp2_types.F:345
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:122
pw_pool_types::pw_pool_type
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:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_dispersion_types::qs_dispersion_type
Definition qs_dispersion_types.F:33
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:217
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:136
qs_linres_types::linres_control_type
General settings for linear response calculations.
Definition qs_linres_types.F:53
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_p_env_types::qs_p_env_type
Represent a qs system that is perturbed. Can calculate the linear operator and the rhs of the system ...
Definition qs_p_env_types.F:58
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
task_list_types::task_list_type
Definition task_list_types.F:58
virial_types::virial_type
Definition virial_types.F:26