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