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