(git:f2099e5)
Loading...
Searching...
No Matches
qs_ks_utils.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief routines that build the Kohn-Sham matrix (i.e calculate the coulomb
10!> and xc parts
11!> \par History
12!> 05.2002 moved from qs_scf (see there the history) [fawzi]
13!> JGH [30.08.02] multi-grid arrays independent from density and potential
14!> 10.2002 introduced pools, uses updated rho as input,
15!> removed most temporary variables, renamed may vars,
16!> began conversion to LSD [fawzi]
17!> 10.2004 moved calculate_w_matrix here [Joost VandeVondele]
18!> introduced energy derivative wrt MOs [Joost VandeVondele]
19!> \author Fawzi Mohamed
20! **************************************************************************************************
21
22MODULE qs_ks_utils
23 USE admm_types, ONLY: admm_type,&
24 get_admm_env
25 USE atomic_kind_types, ONLY: atomic_kind_type
26 USE cell_types, ONLY: cell_type
27 USE cp_control_types, ONLY: dft_control_type
28 USE cp_dbcsr_api, ONLY: &
29 dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_info, dbcsr_init_p, &
30 dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_scale, dbcsr_set, dbcsr_type
31 USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
32 dbcsr_scale_by_vector
33 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
34 copy_fm_to_dbcsr,&
35 cp_dbcsr_plus_fm_fm_t,&
36 cp_dbcsr_sm_fm_multiply,&
37 dbcsr_allocate_matrix_set,&
38 dbcsr_deallocate_matrix_set
39 USE cp_ddapc, ONLY: cp_ddapc_apply_cd
40 USE cp_fm_struct, ONLY: cp_fm_struct_create,&
41 cp_fm_struct_release,&
42 cp_fm_struct_type
43 USE cp_fm_types, ONLY: cp_fm_create,&
44 cp_fm_get_info,&
45 cp_fm_release,&
46 cp_fm_set_all,&
47 cp_fm_to_fm,&
48 cp_fm_type
49 USE cp_log_handling, ONLY: cp_get_default_logger,&
50 cp_logger_type,&
51 cp_to_string
52 USE cp_output_handling, ONLY: cp_p_file,&
53 cp_print_key_finished_output,&
54 cp_print_key_should_output,&
55 cp_print_key_unit_nr
56 USE hfx_admm_utils, ONLY: tddft_hfx_matrix
57 USE hfx_derivatives, ONLY: derivatives_four_center
58 USE hfx_types, ONLY: hfx_type
59 USE input_constants, ONLY: &
60 cdft_alpha_constraint, cdft_beta_constraint, cdft_charge_constraint, &
61 cdft_magnetization_constraint, do_admm_aux_exch_func_none, do_ppl_grid, sic_ad, sic_eo, &
62 sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none
63 USE input_section_types, ONLY: section_vals_get_subs_vals,&
64 section_vals_type,&
65 section_vals_val_get
66 USE kahan_sum, ONLY: accurate_dot_product,&
67 accurate_sum
68 USE kinds, ONLY: default_string_length,&
69 dp
70 USE kpoint_types, ONLY: get_kpoint_info,&
71 kpoint_type
72 USE lri_environment_methods, ONLY: v_int_ppl_update
73 USE lri_environment_types, ONLY: lri_density_type,&
74 lri_environment_type,&
75 lri_kind_type
76 USE lri_forces, ONLY: calculate_lri_forces,&
77 calculate_ri_forces
78 USE lri_ks_methods, ONLY: calculate_lri_ks_matrix,&
79 calculate_ri_ks_matrix
80 USE mathlib, ONLY: abnormal_value
81 USE message_passing, ONLY: mp_para_env_type
82 USE ps_implicit_types, ONLY: mixed_bc,&
83 mixed_periodic_bc,&
84 neumann_bc,&
85 periodic_bc
86 USE pw_env_types, ONLY: pw_env_get,&
87 pw_env_type
88 USE pw_methods, ONLY: pw_axpy,&
89 pw_copy,&
90 pw_integral_ab,&
91 pw_integrate_function,&
92 pw_scale,&
93 pw_transfer,&
94 pw_zero
95 USE pw_poisson_methods, ONLY: pw_poisson_solve
96 USE pw_poisson_types, ONLY: pw_poisson_implicit,&
97 pw_poisson_type
98 USE pw_pool_types, ONLY: pw_pool_type
99 USE pw_types, ONLY: pw_c1d_gs_type,&
100 pw_r3d_rs_type
101 USE qs_cdft_types, ONLY: cdft_control_type
102 USE qs_charges_types, ONLY: qs_charges_type
103 USE qs_collocate_density, ONLY: calculate_rho_elec
104 USE qs_energy_types, ONLY: qs_energy_type
105 USE qs_environment_types, ONLY: get_qs_env,&
106 qs_environment_type
107 USE qs_force_types, ONLY: qs_force_type
108 USE qs_integrate_potential, ONLY: integrate_v_rspace,&
109 integrate_v_rspace_diagonal,&
110 integrate_v_rspace_one_center
111 USE qs_kind_types, ONLY: get_qs_kind_set,&
112 qs_kind_type
113 USE qs_ks_qmmm_methods, ONLY: qmmm_modify_hartree_pot
114 USE qs_ks_types, ONLY: get_ks_env,&
115 qs_ks_env_type
116 USE qs_mo_types, ONLY: get_mo_set,&
117 mo_set_type
118 USE qs_rho_types, ONLY: qs_rho_get,&
119 qs_rho_type
120 USE skala_gpw_functional, ONLY: native_skala_gapw_composite_direct_ao,&
121 native_skala_gapw_composite_reference
122 USE task_list_types, ONLY: task_list_type
123 USE virial_types, ONLY: virial_type
124 USE xc, ONLY: xc_exc_calc,&
125 xc_vxc_pw_create
126#include "./base/base_uses.f90"
127
128 IMPLICIT NONE
129
130 PRIVATE
131
132 LOGICAL, PARAMETER :: debug_this_module = .true.
133 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_utils'
134
135 PUBLIC :: low_spin_roks, sic_explicit_orbitals, calc_v_sic_rspace, check_sum_energies, print_densities, &
136 print_detailed_energy, compute_matrix_vxc, compute_matrix_vxc_kp, sum_up_and_integrate, &
137 calculate_zmp_potential, get_embed_potential_energy
138
139CONTAINS
140
141! **************************************************************************************************
142!> \brief do ROKS calculations yielding low spin states
143!> \param energy ...
144!> \param qs_env ...
145!> \param dft_control ...
146!> \param do_hfx ...
147!> \param just_energy ...
148!> \param calculate_forces ...
149!> \param auxbas_pw_pool ...
150! **************************************************************************************************
151 SUBROUTINE low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
152 calculate_forces, auxbas_pw_pool)
153
154 TYPE(qs_energy_type), POINTER :: energy
155 TYPE(qs_environment_type), POINTER :: qs_env
156 TYPE(dft_control_type), POINTER :: dft_control
157 LOGICAL, INTENT(IN) :: do_hfx, just_energy, calculate_forces
158 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
159
160 CHARACTER(*), PARAMETER :: routinen = 'low_spin_roks'
161
162 INTEGER :: handle, irep, ispin, iterm, k, k_alpha, &
163 k_beta, n_rep, nelectron, nspin, nterms
164 INTEGER, DIMENSION(:), POINTER :: ivec
165 INTEGER, DIMENSION(:, :, :), POINTER :: occupations
166 LOGICAL :: compute_virial, in_range, &
167 uniform_occupation
168 REAL(kind=dp) :: ehfx, exc
169 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
170 REAL(kind=dp), DIMENSION(:), POINTER :: energy_scaling, rvec, scaling
171 TYPE(cell_type), POINTER :: cell
172 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_hfx, matrix_p, mdummy, &
173 mo_derivs, rho_ao
174 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p2
175 TYPE(dbcsr_type), POINTER :: dbcsr_deriv, fm_deriv, fm_scaled, &
176 mo_coeff
177 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
178 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
179 TYPE(mp_para_env_type), POINTER :: para_env
180 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
181 TYPE(pw_env_type), POINTER :: pw_env
182 TYPE(pw_pool_type), POINTER :: xc_pw_pool
183 TYPE(pw_r3d_rs_type) :: work_v_rspace
184 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc, vxc_tau
185 TYPE(pw_r3d_rs_type), POINTER :: weights
186 TYPE(qs_ks_env_type), POINTER :: ks_env
187 TYPE(qs_rho_type), POINTER :: rho
188 TYPE(section_vals_type), POINTER :: hfx_section, input, &
189 low_spin_roks_section, xc_section
190 TYPE(virial_type), POINTER :: virial
191
192 IF (.NOT. dft_control%low_spin_roks) RETURN
193
194 CALL timeset(routinen, handle)
195
196 NULLIFY (ks_env, rho_ao)
197
198 ! Test for not compatible options
199 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
200 CALL cp_abort(__location__, "GAPW/GAPW_XC are not compatible with low spin ROKS method.")
201 END IF
202 IF (dft_control%do_admm) THEN
203 CALL cp_abort(__location__, "ADMM not compatible with low spin ROKS method.")
204 END IF
205 IF (dft_control%do_admm) THEN
206 IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
207 CALL cp_abort(__location__, "ADMM with XC correction functional "// &
208 "not compatible with low spin ROKS method.")
209 END IF
210 END IF
211 IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
212 dft_control%qs_control%xtb) THEN
213 CALL cp_abort(__location__, "SE/xTB/DFTB are not compatible with low spin ROKS method.")
214 END IF
215
216 CALL get_qs_env(qs_env, &
217 ks_env=ks_env, &
218 mo_derivs=mo_derivs, &
219 mos=mo_array, &
220 rho=rho, &
221 pw_env=pw_env, &
222 xcint_weights=weights, &
223 input=input, &
224 cell=cell, &
225 virial=virial)
226
227 CALL qs_rho_get(rho, rho_ao=rho_ao)
228
229 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
230 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
231 hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
232
233 ! No accurate integration possible (as there is no GAPW)
234 IF (ASSOCIATED(weights)) THEN
235 CALL cp_abort(__location__, "No accurate xc integration possible.")
236 END IF
237 ! some assumptions need to be checked
238 ! we have two spins
239 cpassert(SIZE(mo_array, 1) == 2)
240 nspin = 2
241 ! we want uniform occupations
242 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
243 cpassert(uniform_occupation)
244 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation)
245 cpassert(uniform_occupation)
246 IF (do_hfx .AND. calculate_forces .AND. compute_virial) THEN
247 CALL cp_abort(__location__, "ROKS virial with HFX not available.")
248 END IF
249
250 NULLIFY (dbcsr_deriv)
251 CALL dbcsr_init_p(dbcsr_deriv)
252 CALL dbcsr_copy(dbcsr_deriv, mo_derivs(1)%matrix)
253 CALL dbcsr_set(dbcsr_deriv, 0.0_dp)
254
255 ! basic info
256 CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
257 CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_alpha)
258 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff)
259 CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_beta)
260
261 ! read the input
262 low_spin_roks_section => section_vals_get_subs_vals(input, "DFT%LOW_SPIN_ROKS")
263
264 CALL section_vals_val_get(low_spin_roks_section, "ENERGY_SCALING", r_vals=rvec)
265 nterms = SIZE(rvec)
266 ALLOCATE (energy_scaling(nterms))
267 energy_scaling = rvec !? just wondering, should this add up to 1, in which case we should cpp?
268
269 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", n_rep_val=n_rep)
270 cpassert(n_rep == nterms)
271 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=1, i_vals=ivec)
272 nelectron = SIZE(ivec)
273 cpassert(nelectron == k_alpha - k_beta)
274 ALLOCATE (occupations(2, nelectron, nterms))
275 occupations = 0
276 DO iterm = 1, nterms
277 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=iterm, i_vals=ivec)
278 cpassert(nelectron == SIZE(ivec))
279 in_range = all(ivec >= 1) .AND. all(ivec <= 2)
280 cpassert(in_range)
281 DO k = 1, nelectron
282 occupations(ivec(k), k, iterm) = 1
283 END DO
284 END DO
285
286 ! set up general data structures
287 ! density matrices, kohn-sham matrices
288
289 NULLIFY (matrix_p)
290 CALL dbcsr_allocate_matrix_set(matrix_p, nspin)
291 DO ispin = 1, nspin
292 ALLOCATE (matrix_p(ispin)%matrix)
293 CALL dbcsr_copy(matrix_p(ispin)%matrix, rho_ao(1)%matrix, &
294 name="density matrix low spin roks")
295 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
296 END DO
297
298 NULLIFY (matrix_h)
299 CALL dbcsr_allocate_matrix_set(matrix_h, nspin)
300 DO ispin = 1, nspin
301 ALLOCATE (matrix_h(ispin)%matrix)
302 CALL dbcsr_copy(matrix_h(ispin)%matrix, rho_ao(1)%matrix, &
303 name="KS matrix low spin roks")
304 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
305 END DO
306
307 IF (do_hfx) THEN
308 NULLIFY (matrix_hfx)
309 CALL dbcsr_allocate_matrix_set(matrix_hfx, nspin)
310 DO ispin = 1, nspin
311 ALLOCATE (matrix_hfx(ispin)%matrix)
312 CALL dbcsr_copy(matrix_hfx(ispin)%matrix, rho_ao(1)%matrix, &
313 name="HFX matrix low spin roks")
314 END DO
315 END IF
316
317 ! grids in real and g space for rho and vxc
318 ! tau functionals are not supported
319 NULLIFY (tau, vxc_tau, vxc)
320 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
321
322 ALLOCATE (rho_r(nspin))
323 ALLOCATE (rho_g(nspin))
324 DO ispin = 1, nspin
325 CALL auxbas_pw_pool%create_pw(rho_r(ispin))
326 CALL auxbas_pw_pool%create_pw(rho_g(ispin))
327 END DO
328 CALL auxbas_pw_pool%create_pw(work_v_rspace)
329
330 ! get mo matrices needed to construct the density matrices
331 ! we will base all on the alpha spin matrix, obviously possible in ROKS
332 CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
333 NULLIFY (fm_scaled, fm_deriv)
334 CALL dbcsr_init_p(fm_scaled)
335 CALL dbcsr_init_p(fm_deriv)
336 CALL dbcsr_copy(fm_scaled, mo_coeff)
337 CALL dbcsr_copy(fm_deriv, mo_coeff)
338
339 ALLOCATE (scaling(k_alpha))
340
341 ! for each term, add it with the given scaling factor to the energy, and compute the required derivatives
342 DO iterm = 1, nterms
343
344 DO ispin = 1, nspin
345 ! compute the proper density matrices with the required occupations
346 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
347 scaling = 1.0_dp
348 scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
349 CALL dbcsr_copy(fm_scaled, mo_coeff)
350 CALL dbcsr_scale_by_vector(fm_scaled, scaling, side='right')
351 CALL dbcsr_multiply('n', 't', 1.0_dp, mo_coeff, fm_scaled, &
352 0.0_dp, matrix_p(ispin)%matrix, retain_sparsity=.true.)
353 ! compute the densities on the grid
354 CALL calculate_rho_elec(matrix_p=matrix_p(ispin)%matrix, &
355 rho=rho_r(ispin), rho_gspace=rho_g(ispin), &
356 ks_env=ks_env)
357 END DO
358
359 ! compute the exchange energies / potential if needed
360 IF (just_energy) THEN
361 exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
362 weights=weights, pw_pool=xc_pw_pool)
363 ELSE
364 cpassert(.NOT. compute_virial)
365 CALL xc_vxc_pw_create(vxc_rho=vxc, rho_r=rho_r, &
366 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
367 weights=weights, pw_pool=xc_pw_pool, &
368 compute_virial=.false., virial_xc=virial_xc_tmp)
369 END IF
370
371 energy%exc = energy%exc + energy_scaling(iterm)*exc
372
373 IF (do_hfx) THEN
374 ! Add Hartree-Fock contribution
375 DO ispin = 1, nspin
376 CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
377 END DO
378 ehfx = energy%ex
379 CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, &
380 recalc_integrals=.false., update_energy=.true.)
381 energy%ex = ehfx + energy_scaling(iterm)*energy%ex
382 END IF
383
384 ! add the corresponding derivatives to the MO derivatives
385 IF (.NOT. just_energy) THEN
386 ! get the potential in matrix form
387 DO ispin = 1, nspin
388 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
389 ! use a work_v_rspace
390 CALL pw_axpy(vxc(ispin), work_v_rspace, energy_scaling(iterm)*vxc(ispin)%pw_grid%dvol, 0.0_dp)
391 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), &
392 qs_env=qs_env, calculate_forces=calculate_forces)
393 CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
394 END DO
395 DEALLOCATE (vxc)
396
397 IF (do_hfx) THEN
398 ! add HFX contribution
399 DO ispin = 1, nspin
400 CALL dbcsr_add(matrix_h(ispin)%matrix, matrix_hfx(ispin)%matrix, &
401 1.0_dp, energy_scaling(iterm))
402 END DO
403 IF (calculate_forces) THEN
404 CALL get_qs_env(qs_env, x_data=x_data, para_env=para_env)
405 IF (x_data(1, 1)%n_rep_hf /= 1) THEN
406 CALL cp_abort(__location__, "Multiple HFX section forces not compatible "// &
407 "with low spin ROKS method.")
408 END IF
409 IF (x_data(1, 1)%do_hfx_ri) THEN
410 CALL cp_abort(__location__, "HFX_RI forces not compatible with low spin ROKS method.")
411 ELSE
412 irep = 1
413 NULLIFY (mdummy)
414 matrix_p2(1:nspin, 1:1) => matrix_p(1:nspin)
415 CALL derivatives_four_center(qs_env, matrix_p2, mdummy, hfx_section, para_env, &
416 irep, compute_virial, &
417 adiabatic_rescale_factor=energy_scaling(iterm))
418 END IF
419 END IF
420
421 END IF
422
423 ! add this to the mo_derivs, again based on the alpha mo_coeff
424 DO ispin = 1, nspin
425 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, &
426 0.0_dp, dbcsr_deriv, last_column=k_alpha)
427
428 scaling = 1.0_dp
429 scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
430 CALL dbcsr_scale_by_vector(dbcsr_deriv, scaling, side='right')
431 CALL dbcsr_add(mo_derivs(1)%matrix, dbcsr_deriv, 1.0_dp, 1.0_dp)
432 END DO
433
434 END IF
435
436 END DO
437
438 ! release allocated memory
439 DO ispin = 1, nspin
440 CALL auxbas_pw_pool%give_back_pw(rho_r(ispin))
441 CALL auxbas_pw_pool%give_back_pw(rho_g(ispin))
442 END DO
443 DEALLOCATE (rho_r, rho_g)
444 CALL dbcsr_deallocate_matrix_set(matrix_p)
445 CALL dbcsr_deallocate_matrix_set(matrix_h)
446 IF (do_hfx) THEN
447 CALL dbcsr_deallocate_matrix_set(matrix_hfx)
448 END IF
449
450 CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
451
452 CALL dbcsr_release_p(fm_deriv)
453 CALL dbcsr_release_p(fm_scaled)
454
455 DEALLOCATE (occupations)
456 DEALLOCATE (energy_scaling)
457 DEALLOCATE (scaling)
458
459 CALL dbcsr_release_p(dbcsr_deriv)
460
461 CALL timestop(handle)
462
463 END SUBROUTINE low_spin_roks
464! **************************************************************************************************
465!> \brief do sic calculations on explicit orbitals
466!> \param energy ...
467!> \param qs_env ...
468!> \param dft_control ...
469!> \param poisson_env ...
470!> \param just_energy ...
471!> \param calculate_forces ...
472!> \param auxbas_pw_pool ...
473! **************************************************************************************************
474 SUBROUTINE sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, &
475 calculate_forces, auxbas_pw_pool)
476
477 TYPE(qs_energy_type), POINTER :: energy
478 TYPE(qs_environment_type), POINTER :: qs_env
479 TYPE(dft_control_type), POINTER :: dft_control
480 TYPE(pw_poisson_type), POINTER :: poisson_env
481 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
482 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
483
484 CHARACTER(*), PARAMETER :: routinen = 'sic_explicit_orbitals'
485
486 INTEGER :: handle, i, iorb, k_alpha, k_beta, norb
487 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: sic_orbital_list
488 LOGICAL :: compute_virial, uniform_occupation
489 REAL(kind=dp) :: ener, exc
490 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
491 TYPE(cell_type), POINTER :: cell
492 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
493 TYPE(cp_fm_type) :: matrix_hv, matrix_v
494 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_derivs_local
495 TYPE(cp_fm_type), POINTER :: mo_coeff
496 TYPE(dbcsr_p_type) :: orb_density_matrix_p, orb_h_p
497 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, rho_ao, tmp_dbcsr
498 TYPE(dbcsr_type), POINTER :: orb_density_matrix, orb_h
499 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
500 TYPE(pw_c1d_gs_type) :: work_v_gspace
501 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
502 TYPE(pw_c1d_gs_type), TARGET :: orb_rho_g, tmp_g
503 TYPE(pw_env_type), POINTER :: pw_env
504 TYPE(pw_pool_type), POINTER :: xc_pw_pool
505 TYPE(pw_r3d_rs_type) :: work_v_rspace
506 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc, vxc_tau
507 TYPE(pw_r3d_rs_type), POINTER :: weights
508 TYPE(pw_r3d_rs_type), TARGET :: orb_rho_r, tmp_r
509 TYPE(qs_ks_env_type), POINTER :: ks_env
510 TYPE(qs_rho_type), POINTER :: rho
511 TYPE(section_vals_type), POINTER :: input, xc_section
512 TYPE(virial_type), POINTER :: virial
513
514 IF (dft_control%sic_method_id /= sic_eo) RETURN
515
516 CALL timeset(routinen, handle)
517
518 NULLIFY (tau, vxc_tau, mo_derivs, ks_env, rho_ao)
519
520 ! generate the lists of orbitals that need sic treatment
521 CALL get_qs_env(qs_env, &
522 ks_env=ks_env, &
523 mo_derivs=mo_derivs, &
524 mos=mo_array, &
525 rho=rho, &
526 xcint_weights=weights, &
527 pw_env=pw_env, &
528 input=input, &
529 cell=cell, &
530 virial=virial)
531
532 CALL qs_rho_get(rho, rho_ao=rho_ao)
533
534 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
535 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
536
537 DO i = 1, SIZE(mo_array) !fm->dbcsr
538 IF (mo_array(i)%use_mo_coeff_b) THEN !fm->dbcsr
539 CALL copy_dbcsr_to_fm(mo_array(i)%mo_coeff_b, &
540 mo_array(i)%mo_coeff) !fm->dbcsr
541 END IF !fm->dbcsr
542 END DO !fm->dbcsr
543
544 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
545
546 ! we have two spins
547 cpassert(SIZE(mo_array, 1) == 2)
548 ! we want uniform occupations
549 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
550 cpassert(uniform_occupation)
551 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff, uniform_occupation=uniform_occupation)
552 cpassert(uniform_occupation)
553
554 NULLIFY (tmp_dbcsr)
555 CALL dbcsr_allocate_matrix_set(tmp_dbcsr, SIZE(mo_derivs, 1))
556 DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
557 !
558 NULLIFY (tmp_dbcsr(i)%matrix)
559 CALL dbcsr_init_p(tmp_dbcsr(i)%matrix)
560 CALL dbcsr_copy(tmp_dbcsr(i)%matrix, mo_derivs(i)%matrix)
561 CALL dbcsr_set(tmp_dbcsr(i)%matrix, 0.0_dp)
562 END DO !fm->dbcsr
563
564 k_alpha = 0; k_beta = 0
565 SELECT CASE (dft_control%sic_list_id)
566 CASE (sic_list_all)
567
568 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
569 CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
570
571 IF (SIZE(mo_array, 1) > 1) THEN
572 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
573 CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
574 END IF
575
576 norb = k_alpha + k_beta
577 ALLOCATE (sic_orbital_list(3, norb))
578
579 iorb = 0
580 DO i = 1, k_alpha
581 iorb = iorb + 1
582 sic_orbital_list(1, iorb) = 1
583 sic_orbital_list(2, iorb) = i
584 sic_orbital_list(3, iorb) = 1
585 END DO
586 DO i = 1, k_beta
587 iorb = iorb + 1
588 sic_orbital_list(1, iorb) = 2
589 sic_orbital_list(2, iorb) = i
590 IF (SIZE(mo_derivs, 1) == 1) THEN
591 sic_orbital_list(3, iorb) = 1
592 ELSE
593 sic_orbital_list(3, iorb) = 2
594 END IF
595 END DO
596
597 CASE (sic_list_unpaired)
598 ! we have two spins
599 cpassert(SIZE(mo_array, 1) == 2)
600 ! we have them restricted
601 cpassert(SIZE(mo_derivs, 1) == 1)
602 cpassert(dft_control%restricted)
603
604 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
605 CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
606
607 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
608 CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
609
610 norb = k_alpha - k_beta
611 ALLOCATE (sic_orbital_list(3, norb))
612
613 iorb = 0
614 DO i = k_beta + 1, k_alpha
615 iorb = iorb + 1
616 sic_orbital_list(1, iorb) = 1
617 sic_orbital_list(2, iorb) = i
618 ! we are guaranteed to be restricted
619 sic_orbital_list(3, iorb) = 1
620 END DO
621
622 CASE DEFAULT
623 cpabort("Unknown dft_control%sic_list_id")
624 END SELECT
625
626 ! data needed for each of the orbs
627 CALL auxbas_pw_pool%create_pw(orb_rho_r)
628 CALL auxbas_pw_pool%create_pw(tmp_r)
629 CALL auxbas_pw_pool%create_pw(orb_rho_g)
630 CALL auxbas_pw_pool%create_pw(tmp_g)
631 CALL auxbas_pw_pool%create_pw(work_v_gspace)
632 CALL auxbas_pw_pool%create_pw(work_v_rspace)
633
634 ALLOCATE (orb_density_matrix)
635 CALL dbcsr_copy(orb_density_matrix, rho_ao(1)%matrix, &
636 name="orb_density_matrix")
637 CALL dbcsr_set(orb_density_matrix, 0.0_dp)
638 orb_density_matrix_p%matrix => orb_density_matrix
639
640 ALLOCATE (orb_h)
641 CALL dbcsr_copy(orb_h, rho_ao(1)%matrix, &
642 name="orb_density_matrix")
643 CALL dbcsr_set(orb_h, 0.0_dp)
644 orb_h_p%matrix => orb_h
645
646 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
647
648 CALL cp_fm_struct_create(fm_struct_tmp, ncol_global=1, &
649 template_fmstruct=mo_coeff%matrix_struct)
650 CALL cp_fm_create(matrix_v, fm_struct_tmp, name="matrix_v")
651 CALL cp_fm_create(matrix_hv, fm_struct_tmp, name="matrix_hv")
652 CALL cp_fm_struct_release(fm_struct_tmp)
653
654 ALLOCATE (mo_derivs_local(SIZE(mo_array, 1)))
655 DO i = 1, SIZE(mo_array, 1)
656 CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff)
657 CALL cp_fm_create(mo_derivs_local(i), mo_coeff%matrix_struct)
658 END DO
659
660 ALLOCATE (rho_r(2))
661 rho_r(1) = orb_rho_r
662 rho_r(2) = tmp_r
663 CALL pw_zero(tmp_r)
664
665 ALLOCATE (rho_g(2))
666 rho_g(1) = orb_rho_g
667 rho_g(2) = tmp_g
668 CALL pw_zero(tmp_g)
669
670 NULLIFY (vxc)
671 ! now apply to SIC correction to each selected orbital
672 DO iorb = 1, norb
673 ! extract the proper orbital from the mo_coeff
674 CALL get_mo_set(mo_set=mo_array(sic_orbital_list(1, iorb)), mo_coeff=mo_coeff)
675 CALL cp_fm_to_fm(mo_coeff, matrix_v, 1, sic_orbital_list(2, iorb), 1)
676
677 ! construct the density matrix and the corresponding density
678 CALL dbcsr_set(orb_density_matrix, 0.0_dp)
679 CALL cp_dbcsr_plus_fm_fm_t(orb_density_matrix, matrix_v=matrix_v, ncol=1, &
680 alpha=1.0_dp)
681
682 CALL calculate_rho_elec(matrix_p=orb_density_matrix, &
683 rho=orb_rho_r, rho_gspace=orb_rho_g, &
684 ks_env=ks_env)
685
686 ! compute the energy functional for this orbital and its derivative
687
688 CALL pw_poisson_solve(poisson_env, orb_rho_g, ener, work_v_gspace)
689 ! no PBC correction is done here, see "calc_v_sic_rspace" for SIC methods
690 ! with PBC aware corrections
691 energy%hartree = energy%hartree - dft_control%sic_scaling_a*ener
692 IF (.NOT. just_energy) THEN
693 CALL pw_transfer(work_v_gspace, work_v_rspace)
694 CALL pw_scale(work_v_rspace, -dft_control%sic_scaling_a*work_v_rspace%pw_grid%dvol)
695 CALL dbcsr_set(orb_h, 0.0_dp)
696 END IF
697
698 IF (just_energy) THEN
699 exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
700 weights=weights, pw_pool=xc_pw_pool)
701 ELSE
702 cpassert(.NOT. compute_virial)
703 CALL xc_vxc_pw_create(vxc_rho=vxc, rho_r=rho_r, &
704 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
705 weights=weights, pw_pool=xc_pw_pool, &
706 compute_virial=compute_virial, virial_xc=virial_xc_tmp)
707 ! add to the existing work_v_rspace
708 CALL pw_axpy(vxc(1), work_v_rspace, -dft_control%sic_scaling_b*vxc(1)%pw_grid%dvol)
709 END IF
710 energy%exc = energy%exc - dft_control%sic_scaling_b*exc
711
712 IF (.NOT. just_energy) THEN
713 ! note, orb_h (which is being pointed to with orb_h_p) is zeroed above
714 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=orb_density_matrix_p, hmat=orb_h_p, &
715 qs_env=qs_env, calculate_forces=calculate_forces)
716
717 ! add this to the mo_derivs
718 CALL cp_dbcsr_sm_fm_multiply(orb_h, matrix_v, matrix_hv, 1)
719 ! silly trick, copy to an array of the right size and add to mo_derivs
720 CALL cp_fm_set_all(mo_derivs_local(sic_orbital_list(3, iorb)), 0.0_dp)
721 CALL cp_fm_to_fm(matrix_hv, mo_derivs_local(sic_orbital_list(3, iorb)), 1, 1, sic_orbital_list(2, iorb))
722 CALL copy_fm_to_dbcsr(mo_derivs_local(sic_orbital_list(3, iorb)), &
723 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix)
724 CALL dbcsr_add(mo_derivs(sic_orbital_list(3, iorb))%matrix, &
725 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix, 1.0_dp, 1.0_dp)
726 !
727 ! need to deallocate vxc
728 CALL xc_pw_pool%give_back_pw(vxc(1))
729 CALL xc_pw_pool%give_back_pw(vxc(2))
730 DEALLOCATE (vxc)
731
732 END IF
733
734 END DO
735
736 CALL auxbas_pw_pool%give_back_pw(orb_rho_r)
737 CALL auxbas_pw_pool%give_back_pw(tmp_r)
738 CALL auxbas_pw_pool%give_back_pw(orb_rho_g)
739 CALL auxbas_pw_pool%give_back_pw(tmp_g)
740 CALL auxbas_pw_pool%give_back_pw(work_v_gspace)
741 CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
742
743 CALL dbcsr_deallocate_matrix(orb_density_matrix)
744 CALL dbcsr_deallocate_matrix(orb_h)
745 CALL cp_fm_release(matrix_v)
746 CALL cp_fm_release(matrix_hv)
747 CALL cp_fm_release(mo_derivs_local)
748 DEALLOCATE (rho_r)
749 DEALLOCATE (rho_g)
750
751 CALL dbcsr_deallocate_matrix_set(tmp_dbcsr) !fm->dbcsr
752
753 CALL timestop(handle)
754
755 END SUBROUTINE sic_explicit_orbitals
756
757! **************************************************************************************************
758!> \brief do sic calculations on the spin density
759!> \param v_sic_rspace ...
760!> \param energy ...
761!> \param qs_env ...
762!> \param dft_control ...
763!> \param rho ...
764!> \param poisson_env ...
765!> \param just_energy ...
766!> \param calculate_forces ...
767!> \param auxbas_pw_pool ...
768! **************************************************************************************************
769 SUBROUTINE calc_v_sic_rspace(v_sic_rspace, energy, &
770 qs_env, dft_control, rho, poisson_env, just_energy, &
771 calculate_forces, auxbas_pw_pool)
772
773 TYPE(pw_r3d_rs_type), POINTER :: v_sic_rspace
774 TYPE(qs_energy_type), POINTER :: energy
775 TYPE(qs_environment_type), POINTER :: qs_env
776 TYPE(dft_control_type), POINTER :: dft_control
777 TYPE(qs_rho_type), POINTER :: rho
778 TYPE(pw_poisson_type), POINTER :: poisson_env
779 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
780 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
781
782 INTEGER :: i, nelec, nelec_a, nelec_b, nforce
783 REAL(kind=dp) :: ener, full_scaling, scaling
784 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: store_forces
785 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
786 TYPE(pw_c1d_gs_type) :: work_rho, work_v
787 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
788 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
789
790 NULLIFY (mo_array, rho_g)
791
792 IF (dft_control%sic_method_id == sic_none) RETURN
793 IF (dft_control%sic_method_id == sic_eo) RETURN
794
795 IF (dft_control%qs_control%gapw) THEN
796 cpabort("sic and GAPW not yet compatible")
797 END IF
798
799 ! OK, right now we like two spins to do sic, could be relaxed for AD
800 cpassert(dft_control%nspins == 2)
801
802 CALL auxbas_pw_pool%create_pw(work_rho)
803 CALL auxbas_pw_pool%create_pw(work_v)
804
805 CALL qs_rho_get(rho, rho_g=rho_g)
806
807 ! Hartree sic corrections
808 SELECT CASE (dft_control%sic_method_id)
809 CASE (sic_mauri_us, sic_mauri_spz)
810 CALL pw_copy(rho_g(1), work_rho)
811 CALL pw_axpy(rho_g(2), work_rho, alpha=-1._dp)
812 CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
813 CASE (sic_ad)
814 ! find out how many elecs we have
815 CALL get_qs_env(qs_env, mos=mo_array)
816 CALL get_mo_set(mo_set=mo_array(1), nelectron=nelec_a)
817 CALL get_mo_set(mo_set=mo_array(2), nelectron=nelec_b)
818 nelec = nelec_a + nelec_b
819 CALL pw_copy(rho_g(1), work_rho)
820 CALL pw_axpy(rho_g(2), work_rho)
821 scaling = 1.0_dp/real(nelec, kind=dp)
822 CALL pw_scale(work_rho, scaling)
823 CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
824 CASE DEFAULT
825 cpabort("Unknown sic method id")
826 END SELECT
827
828 ! Correct for DDAP charges (if any)
829 ! storing whatever force might be there from previous decoupling
830 IF (calculate_forces) THEN
831 CALL get_qs_env(qs_env=qs_env, force=force)
832 nforce = 0
833 DO i = 1, SIZE(force)
834 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
835 END DO
836 ALLOCATE (store_forces(3, nforce))
837 nforce = 0
838 DO i = 1, SIZE(force)
839 store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2)) = force(i)%ch_pulay(:, :)
840 force(i)%ch_pulay(:, :) = 0.0_dp
841 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
842 END DO
843 END IF
844
845 CALL cp_ddapc_apply_cd(qs_env, &
846 work_rho, &
847 ener, &
848 v_hartree_gspace=work_v, &
849 calculate_forces=calculate_forces, &
850 itype_of_density="SPIN")
851
852 SELECT CASE (dft_control%sic_method_id)
853 CASE (sic_mauri_us, sic_mauri_spz)
854 full_scaling = -dft_control%sic_scaling_a
855 CASE (sic_ad)
856 full_scaling = -dft_control%sic_scaling_a*nelec
857 CASE DEFAULT
858 cpabort("Unknown sic method id")
859 END SELECT
860 energy%hartree = energy%hartree + full_scaling*ener
861
862 ! add scaled forces, restoring the old
863 IF (calculate_forces) THEN
864 nforce = 0
865 DO i = 1, SIZE(force)
866 force(i)%ch_pulay(:, :) = force(i)%ch_pulay(:, :)*full_scaling + &
867 store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2))
868 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
869 END DO
870 END IF
871
872 IF (.NOT. just_energy) THEN
873 ALLOCATE (v_sic_rspace)
874 CALL auxbas_pw_pool%create_pw(v_sic_rspace)
875 CALL pw_transfer(work_v, v_sic_rspace)
876 ! also take into account the scaling (in addition to the volume element)
877 CALL pw_scale(v_sic_rspace, &
878 dft_control%sic_scaling_a*v_sic_rspace%pw_grid%dvol)
879 END IF
880
881 CALL auxbas_pw_pool%give_back_pw(work_rho)
882 CALL auxbas_pw_pool%give_back_pw(work_v)
883
884 END SUBROUTINE calc_v_sic_rspace
885
886! **************************************************************************************************
887!> \brief Check each term of the energy and sum to total
888!> \param energy ...
889! **************************************************************************************************
890 SUBROUTINE check_sum_energies(energy)
891 TYPE(qs_energy_type), INTENT(INOUT), POINTER :: energy
892
893 ! Reset the total energy counter, just in case
894 ! it is used to store some intermediate value
895 energy%total = 0.0_dp
896
897 CALL check_sum_energy(energy%total, energy%core_overlap, "Core overlap energy")
898 CALL check_sum_energy(energy%total, energy%core_self, "Core self energy")
899 CALL check_sum_energy(energy%total, energy%core_cneo, "Quantum nuclear core energy")
900 CALL check_sum_energy(energy%total, energy%core, "Core Hamiltonian energy")
901 CALL check_sum_energy(energy%total, energy%hartree, "Hartree energy")
902 CALL check_sum_energy(energy%total, energy%hartree_1c, "Energy from GAPW local Eh = 1 center integrals")
903 CALL check_sum_energy(energy%total, energy%exc, "Exchange-correlation energy")
904 CALL check_sum_energy(energy%total, energy%exc1, "Exc1 energy")
905 CALL check_sum_energy(energy%total, energy%ex, "Ex energy")
906 CALL check_sum_energy(energy%total, energy%dispersion, "Dispersion energy")
907 CALL check_sum_energy(energy%total, energy%gcp, "gCP energy")
908 CALL check_sum_energy(energy%total, energy%qmmm_el, "QM/MM Electrostatic energy")
909 CALL check_sum_energy(energy%total, energy%mulliken, "Mulliken restraint energy")
910 CALL check_sum_energy(energy%total, sum(energy%ddapc_restraint), "DDAPC restraint energy")
911 CALL check_sum_energy(energy%total, energy%s2_restraint, "S2 restraint energy")
912 CALL check_sum_energy(energy%total, energy%dft_plus_u, "DFT+U energy")
913 CALL check_sum_energy(energy%total, energy%kTS, "Electronic entropic contribution energy")
914 CALL check_sum_energy(energy%total, energy%efield, "Electric field interaction energy")
915 CALL check_sum_energy(energy%total, energy%efield_core, "Electric field core interaction energy")
916 CALL check_sum_energy(energy%total, energy%ee, "External potential interaction energy")
917 CALL check_sum_energy(energy%total, energy%ee_core, "External potential core interaction energy")
918 CALL check_sum_energy(energy%total, energy%exc_aux_fit, "Wfn fit exchange-correlation energy")
919 CALL check_sum_energy(energy%total, energy%image_charge, "QM/MM image charge energy energy")
920 CALL check_sum_energy(energy%total, energy%sccs_pol, "SCCS polarisation energy")
921 CALL check_sum_energy(energy%total, energy%cdft, "CDFT constraint energy")
922 CALL check_sum_energy(energy%total, energy%exc1_aux_fit, "Wfn fit soft/hard atomic rho1 Exc contribution energy")
923 CALL check_sum_energy(energy%total, energy%embed_corr, "Embedding potential correction energy")
924
925 IF (abnormal_value(energy%total)) THEN
926 CALL cp_abort(__location__, &
927 "Total energy after summing up the terms is an abnormal value (NaN/Inf).")
928 END IF
929
930 END SUBROUTINE check_sum_energies
931
932! **************************************************************************************************
933!> \brief Check one term of the energy and sum to total if normal
934!> \param total ...
935!> \param term ...
936!> \param name ...
937! **************************************************************************************************
938 SUBROUTINE check_sum_energy(total, term, name)
939 REAL(kind=dp), INTENT(INOUT) :: total
940 REAL(kind=dp), INTENT(IN) :: term
941 CHARACTER(LEN=*), INTENT(IN) :: name
942
943 IF (abnormal_value(term)) THEN
944 CALL cp_abort(__location__, &
945 trim(name)//" is an abnormal value (NaN/Inf).")
946 ELSE
947 total = total + term
948 END IF
949
950 END SUBROUTINE check_sum_energy
951
952! **************************************************************************************************
953!> \brief ...
954!> \param qs_env ...
955!> \param rho ...
956! **************************************************************************************************
957 SUBROUTINE print_densities(qs_env, rho)
958 TYPE(qs_environment_type), POINTER :: qs_env
959 TYPE(qs_rho_type), POINTER :: rho
960
961 INTEGER :: img, ispin, n_electrons, output_unit
962 REAL(dp) :: tot1_h, tot1_s, tot_rho_r, trace, &
963 trace_tmp
964 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r_arr
965 TYPE(cell_type), POINTER :: cell
966 TYPE(cp_logger_type), POINTER :: logger
967 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, rho_ao
968 TYPE(dft_control_type), POINTER :: dft_control
969 TYPE(qs_charges_type), POINTER :: qs_charges
970 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
971 TYPE(section_vals_type), POINTER :: input, scf_section
972
973 NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, &
974 dft_control, tot_rho_r_arr, rho_ao)
975
976 logger => cp_get_default_logger()
977
978 CALL get_qs_env(qs_env, &
979 qs_kind_set=qs_kind_set, &
980 cell=cell, qs_charges=qs_charges, &
981 input=input, &
982 matrix_s_kp=matrix_s, &
983 dft_control=dft_control)
984
985 CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
986
987 scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
988 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", &
989 extension=".scfLog")
990
991 CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao)
992 n_electrons = n_electrons - dft_control%charge
993 tot_rho_r = accurate_sum(tot_rho_r_arr)
994
995 trace = 0
996 IF (btest(cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES"), cp_p_file)) THEN
997 DO ispin = 1, dft_control%nspins
998 DO img = 1, dft_control%nimages
999 CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp)
1000 trace = trace + trace_tmp
1001 END DO
1002 END DO
1003 END IF
1004
1005 IF (output_unit > 0) THEN
1006 WRITE (unit=output_unit, fmt="(/,T3,A,T41,F20.10)") "Trace(PS):", trace
1007 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
1008 "Electronic density on regular grids: ", &
1009 tot_rho_r, &
1010 tot_rho_r + &
1011 REAL(n_electrons, dp), &
1012 "Core density on regular grids:", &
1013 qs_charges%total_rho_core_rspace, &
1014 qs_charges%total_rho_core_rspace + &
1015 qs_charges%total_rho1_hard_nuc - &
1016 REAL(n_electrons + dft_control%charge, dp)
1017 END IF
1018 IF (dft_control%qs_control%gapw) THEN
1019 tot1_h = qs_charges%total_rho1_hard(1)
1020 tot1_s = qs_charges%total_rho1_soft(1)
1021 DO ispin = 2, dft_control%nspins
1022 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
1023 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
1024 END DO
1025 IF (output_unit > 0) THEN
1026 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
1027 "Hard and soft densities (Lebedev):", &
1028 tot1_h, tot1_s
1029 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
1030 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
1031 tot_rho_r + tot1_h - tot1_s, &
1032 "Total charge density (r-space): ", &
1033 tot_rho_r + tot1_h - tot1_s &
1034 + qs_charges%total_rho_core_rspace &
1035 + qs_charges%total_rho1_hard_nuc
1036 IF (qs_charges%total_rho1_hard_nuc /= 0.0_dp) THEN
1037 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
1038 "Total CNEO nuc. char. den. (Lebedev): ", &
1039 qs_charges%total_rho1_hard_nuc, &
1040 "Total CNEO soft char. den. (Lebedev): ", &
1041 qs_charges%total_rho1_soft_nuc_lebedev, &
1042 "Total CNEO soft char. den. (r-space): ", &
1043 qs_charges%total_rho1_soft_nuc_rspace, &
1044 "Total soft Rho_e+n+0 (g-space):", &
1045 qs_charges%total_rho_gspace
1046 ELSE
1047 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
1048 "Total Rho_soft + Rho0_soft (g-space):", &
1049 qs_charges%total_rho_gspace
1050 END IF
1051 END IF
1052 qs_charges%background = tot_rho_r + tot1_h - tot1_s + &
1053 qs_charges%total_rho_core_rspace + &
1054 qs_charges%total_rho1_hard_nuc
1055 ! only add total_rho1_hard_nuc for gapw as cneo requires gapw
1056 ELSE IF (dft_control%qs_control%gapw_xc) THEN
1057 tot1_h = qs_charges%total_rho1_hard(1)
1058 tot1_s = qs_charges%total_rho1_soft(1)
1059 DO ispin = 2, dft_control%nspins
1060 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
1061 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
1062 END DO
1063 IF (output_unit > 0) THEN
1064 WRITE (unit=output_unit, fmt="(/,(T3,A,T41,2F20.10))") &
1065 "Hard and soft densities (Lebedev):", &
1066 tot1_h, tot1_s
1067 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
1068 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
1069 accurate_sum(tot_rho_r_arr) + tot1_h - tot1_s
1070 END IF
1071 qs_charges%background = tot_rho_r + &
1072 qs_charges%total_rho_core_rspace
1073 ELSE
1074 IF (output_unit > 0) THEN
1075 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
1076 "Total charge density on r-space grids: ", &
1077 tot_rho_r + &
1078 qs_charges%total_rho_core_rspace, &
1079 "Total charge density g-space grids: ", &
1080 qs_charges%total_rho_gspace
1081 END IF
1082 qs_charges%background = tot_rho_r + &
1083 qs_charges%total_rho_core_rspace
1084 END IF
1085 IF (output_unit > 0) WRITE (unit=output_unit, fmt="()")
1086 qs_charges%background = qs_charges%background/cell%deth
1087
1088 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1089 "PRINT%TOTAL_DENSITIES")
1090
1091 END SUBROUTINE print_densities
1092
1093! **************************************************************************************************
1094!> \brief Print detailed energies
1095!>
1096!> \param qs_env ...
1097!> \param dft_control ...
1098!> \param input ...
1099!> \param energy ...
1100!> \param mulliken_order_p ...
1101!> \par History
1102!> refactoring 04.03.2011 [MI]
1103!> \author
1104! **************************************************************************************************
1105 SUBROUTINE print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
1106
1107 TYPE(qs_environment_type), POINTER :: qs_env
1108 TYPE(dft_control_type), POINTER :: dft_control
1109 TYPE(section_vals_type), POINTER :: input
1110 TYPE(qs_energy_type), POINTER :: energy
1111 REAL(kind=dp), INTENT(IN) :: mulliken_order_p
1112
1113 INTEGER :: bc, n, output_unit, psolver
1114 REAL(kind=dp) :: ddapc_order_p, implicit_ps_ehartree, &
1115 s2_order_p
1116 TYPE(cp_logger_type), POINTER :: logger
1117 TYPE(pw_env_type), POINTER :: pw_env
1118
1119 logger => cp_get_default_logger()
1120
1121 NULLIFY (pw_env)
1122 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1123 psolver = pw_env%poisson_env%parameters%solver
1124
1125 output_unit = cp_print_key_unit_nr(logger, input, "DFT%SCF%PRINT%DETAILED_ENERGY", &
1126 extension=".scfLog")
1127 IF (output_unit > 0) THEN
1128 IF (dft_control%do_admm) THEN
1129 WRITE (unit=output_unit, fmt="((T3,A,T60,F20.10))") &
1130 "Wfn fit exchange-correlation energy: ", energy%exc_aux_fit
1131 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1132 WRITE (unit=output_unit, fmt="((T3,A,T60,F20.10))") &
1133 "Wfn fit soft/hard atomic rho1 Exc contribution: ", energy%exc1_aux_fit
1134 END IF
1135 END IF
1136 IF (dft_control%do_admm) THEN
1137 IF (psolver == pw_poisson_implicit) THEN
1138 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1139 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1140 SELECT CASE (bc)
1141 CASE (mixed_periodic_bc, mixed_bc)
1142 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1143 "Core Hamiltonian energy: ", energy%core, &
1144 "Hartree energy: ", implicit_ps_ehartree, &
1145 "Electric enthalpy: ", energy%hartree, &
1146 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1147 CASE (periodic_bc, neumann_bc)
1148 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1149 "Core Hamiltonian energy: ", energy%core, &
1150 "Hartree energy: ", energy%hartree, &
1151 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1152 END SELECT
1153 ELSE
1154 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1155 "Core Hamiltonian energy: ", energy%core, &
1156 "Hartree energy: ", energy%hartree, &
1157 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1158 END IF
1159 ELSE
1160 !ZMP to print some variables at each step
1161 IF (dft_control%apply_external_density) THEN
1162 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1163 "DOING ZMP CALCULATION FROM EXTERNAL DENSITY "
1164 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1165 "Core Hamiltonian energy: ", energy%core, &
1166 "Hartree energy: ", energy%hartree
1167 ELSE IF (dft_control%apply_external_vxc) THEN
1168 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1169 "DOING ZMP READING EXTERNAL VXC "
1170 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1171 "Core Hamiltonian energy: ", energy%core, &
1172 "Hartree energy: ", energy%hartree
1173 ELSE
1174 IF (psolver == pw_poisson_implicit) THEN
1175 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1176 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1177 SELECT CASE (bc)
1178 CASE (mixed_periodic_bc, mixed_bc)
1179 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1180 "Core Hamiltonian energy: ", energy%core, &
1181 "Hartree energy: ", implicit_ps_ehartree, &
1182 "Electric enthalpy: ", energy%hartree, &
1183 "Exchange-correlation energy: ", energy%exc
1184 CASE (periodic_bc, neumann_bc)
1185 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1186 "Core Hamiltonian energy: ", energy%core, &
1187 "Hartree energy: ", energy%hartree, &
1188 "Exchange-correlation energy: ", energy%exc
1189 END SELECT
1190 ELSE
1191 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1192 "Core Hamiltonian energy: ", energy%core, &
1193 "Hartree energy: ", energy%hartree, &
1194 "Exchange-correlation energy: ", energy%exc
1195 END IF
1196 END IF
1197 END IF
1198
1199 IF (dft_control%apply_external_density) THEN
1200 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1201 "Integral of the (density * v_xc): ", energy%exc
1202 END IF
1203
1204 IF (energy%e_hartree /= 0.0_dp) THEN
1205 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1206 "Coulomb (electron-electron) energy: ", energy%e_hartree
1207 END IF
1208 IF (energy%dispersion /= 0.0_dp) THEN
1209 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1210 "Dispersion energy: ", energy%dispersion
1211 END IF
1212 IF (energy%efield /= 0.0_dp) THEN
1213 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1214 "Electric field interaction energy: ", energy%efield
1215 END IF
1216 IF (energy%gcp /= 0.0_dp) THEN
1217 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1218 "gCP energy: ", energy%gcp
1219 END IF
1220 IF (dft_control%qs_control%gapw) THEN
1221 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1222 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit, &
1223 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
1224 END IF
1225 IF (dft_control%qs_control%gapw_xc) THEN
1226 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1227 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit
1228 END IF
1229 IF (dft_control%dft_plus_u) THEN
1230 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1231 "DFT+U energy:", energy%dft_plus_u
1232 END IF
1233 IF (qs_env%qmmm) THEN
1234 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1235 "QM/MM Electrostatic energy: ", energy%qmmm_el
1236 IF (qs_env%qmmm_env_qm%image_charge) THEN
1237 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1238 "QM/MM image charge energy: ", energy%image_charge
1239 END IF
1240 END IF
1241 IF (dft_control%qs_control%mulliken_restraint) THEN
1242 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1243 "Mulliken restraint (order_p,energy) : ", mulliken_order_p, energy%mulliken
1244 END IF
1245 IF (dft_control%qs_control%ddapc_restraint) THEN
1246 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
1247 ddapc_order_p = &
1248 dft_control%qs_control%ddapc_restraint_control(n)%ddapc_order_p
1249 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1250 "DDAPC restraint (order_p,energy) : ", ddapc_order_p, energy%ddapc_restraint(n)
1251 END DO
1252 END IF
1253 IF (dft_control%qs_control%s2_restraint) THEN
1254 s2_order_p = dft_control%qs_control%s2_restraint_control%s2_order_p
1255 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1256 "S2 restraint (order_p,energy) : ", s2_order_p, energy%s2_restraint
1257 END IF
1258 IF (energy%core_cneo /= 0.0_dp) THEN
1259 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1260 "CNEO| quantum nuclear core energy: ", energy%core_cneo
1261 END IF
1262
1263 END IF ! output_unit
1264 CALL cp_print_key_finished_output(output_unit, logger, input, &
1265 "DFT%SCF%PRINT%DETAILED_ENERGY")
1266
1267 END SUBROUTINE print_detailed_energy
1268
1269! **************************************************************************************************
1270!> \brief compute matrix_vxc, defined via the potential created by qs_vxc_create
1271!> ignores things like tau functional, gapw, sic, ...
1272!> so only OK for GGA & GPW right now
1273!> \param qs_env ...
1274!> \param v_rspace ...
1275!> \param matrix_vxc ...
1276!> \param gapw_full_basis ...
1277!> \par History
1278!> created 23.10.2012 [Joost VandeVondele]
1279!> \author
1280! **************************************************************************************************
1281 SUBROUTINE compute_matrix_vxc(qs_env, v_rspace, matrix_vxc, gapw_full_basis)
1282 TYPE(qs_environment_type), POINTER :: qs_env
1283 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
1284 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
1285 LOGICAL, INTENT(IN), OPTIONAL :: gapw_full_basis
1286
1287 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_matrix_vxc'
1288
1289 INTEGER :: handle, ispin
1290 LOGICAL :: gapw
1291 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1292 TYPE(dft_control_type), POINTER :: dft_control
1293
1294 CALL timeset(routinen, handle)
1295
1296 ! create the matrix using matrix_ks as a template
1297 IF (ASSOCIATED(matrix_vxc)) THEN
1298 CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1299 END IF
1300 CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
1301 ALLOCATE (matrix_vxc(SIZE(matrix_ks)))
1302 DO ispin = 1, SIZE(matrix_ks)
1303 NULLIFY (matrix_vxc(ispin)%matrix)
1304 CALL dbcsr_init_p(matrix_vxc(ispin)%matrix)
1305 CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, &
1306 name="Matrix VXC of spin "//cp_to_string(ispin))
1307 CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp)
1308 END DO
1309
1310 ! and integrate
1311 CALL get_qs_env(qs_env, dft_control=dft_control)
1312 gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1313 IF (PRESENT(gapw_full_basis)) THEN
1314 IF (gapw_full_basis) gapw = .false.
1315 END IF
1316 DO ispin = 1, SIZE(matrix_ks)
1317 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1318 hmat=matrix_vxc(ispin), &
1319 qs_env=qs_env, &
1320 calculate_forces=.false., &
1321 gapw=gapw)
1322 ! scale by the volume element... should really become part of integrate_v_rspace
1323 CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw_grid%dvol)
1324 END DO
1325
1326 CALL timestop(handle)
1327
1328 END SUBROUTINE compute_matrix_vxc
1329
1330! **************************************************************************************************
1331!> \brief Build the XC potential matrix for k-point/image-resolved KS matrices.
1332!> \param qs_env Quickstep environment
1333!> \param v_rspace XC potential on the real-space grid
1334!> \param matrix_vxc_kp k-point/image-resolved XC potential matrix
1335!> \param gapw_full_basis ...
1336!> \author
1337! **************************************************************************************************
1338 SUBROUTINE compute_matrix_vxc_kp(qs_env, v_rspace, matrix_vxc_kp, gapw_full_basis)
1339 TYPE(qs_environment_type), POINTER :: qs_env
1340 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
1341 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_vxc_kp
1342 LOGICAL, INTENT(IN), OPTIONAL :: gapw_full_basis
1343
1344 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_matrix_vxc_kp'
1345
1346 INTEGER :: handle, img, ispin
1347 LOGICAL :: gapw
1348 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat
1349 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
1350 TYPE(dft_control_type), POINTER :: dft_control
1351
1352 CALL timeset(routinen, handle)
1353
1354 ! create the matrix using matrix_ks as a template
1355 IF (ASSOCIATED(matrix_vxc_kp)) THEN
1356 CALL dbcsr_deallocate_matrix_set(matrix_vxc_kp)
1357 END IF
1358 CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks_kp=matrix_ks_kp)
1359 ALLOCATE (matrix_vxc_kp(SIZE(matrix_ks_kp, 1), SIZE(matrix_ks_kp, 2)))
1360 DO img = 1, SIZE(matrix_ks_kp, 2)
1361 DO ispin = 1, SIZE(matrix_ks_kp, 1)
1362 NULLIFY (matrix_vxc_kp(ispin, img)%matrix)
1363 CALL dbcsr_init_p(matrix_vxc_kp(ispin, img)%matrix)
1364 CALL dbcsr_copy(matrix_vxc_kp(ispin, img)%matrix, matrix_ks_kp(ispin, img)%matrix, &
1365 name="Matrix VXC of spin "//cp_to_string(ispin)//" image "//cp_to_string(img))
1366 CALL dbcsr_set(matrix_vxc_kp(ispin, img)%matrix, 0.0_dp)
1367 END DO
1368 END DO
1369
1370 ! and integrate
1371 gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1372 IF (PRESENT(gapw_full_basis)) THEN
1373 IF (gapw_full_basis) gapw = .false.
1374 END IF
1375 DO ispin = 1, SIZE(matrix_ks_kp, 1)
1376 ksmat => matrix_vxc_kp(ispin, :)
1377 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1378 hmat_kp=ksmat, &
1379 qs_env=qs_env, &
1380 calculate_forces=.false., &
1381 gapw=gapw)
1382 ! scale by the volume element... should really become part of integrate_v_rspace
1383 DO img = 1, SIZE(matrix_ks_kp, 2)
1384 CALL dbcsr_scale(matrix_vxc_kp(ispin, img)%matrix, v_rspace(ispin)%pw_grid%dvol)
1385 END DO
1386 END DO
1387
1388 CALL timestop(handle)
1389
1390 END SUBROUTINE compute_matrix_vxc_kp
1391
1392! **************************************************************************************************
1393!> \brief Sum up all potentials defined on the grid and integrate
1394!>
1395!> \param qs_env ...
1396!> \param ks_matrix ...
1397!> \param rho ...
1398!> \param my_rho ...
1399!> \param vppl_rspace ...
1400!> \param v_rspace_new ...
1401!> \param v_rspace_new_aux_fit ...
1402!> \param v_tau_rspace ...
1403!> \param v_tau_rspace_aux_fit ...
1404!> \param v_sic_rspace ...
1405!> \param v_spin_ddapc_rest_r ...
1406!> \param v_sccs_rspace ...
1407!> \param v_rspace_embed ...
1408!> \param cdft_control ...
1409!> \param calculate_forces ...
1410!> \par History
1411!> - refactoring 04.03.2011 [MI]
1412!> - SCCS implementation (16.10.2013,MK)
1413!> \author
1414! **************************************************************************************************
1415 SUBROUTINE sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, &
1416 vppl_rspace, v_rspace_new, &
1417 v_rspace_new_aux_fit, v_tau_rspace, &
1418 v_tau_rspace_aux_fit, &
1419 v_sic_rspace, v_spin_ddapc_rest_r, &
1420 v_sccs_rspace, v_rspace_embed, cdft_control, &
1421 calculate_forces)
1422
1423 TYPE(qs_environment_type), POINTER :: qs_env
1424 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
1425 TYPE(qs_rho_type), POINTER :: rho
1426 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: my_rho
1427 TYPE(pw_r3d_rs_type), POINTER :: vppl_rspace
1428 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_rspace_new_aux_fit, &
1429 v_tau_rspace, v_tau_rspace_aux_fit
1430 TYPE(pw_r3d_rs_type), POINTER :: v_sic_rspace, v_spin_ddapc_rest_r, &
1431 v_sccs_rspace
1432 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_embed
1433 TYPE(cdft_control_type), POINTER :: cdft_control
1434 LOGICAL, INTENT(in) :: calculate_forces
1435
1436 CHARACTER(LEN=*), PARAMETER :: routinen = 'sum_up_and_integrate'
1437
1438 CHARACTER(LEN=default_string_length) :: basis_type
1439 INTEGER :: handle, igroup, ikind, img, ispin, &
1440 nkind, nspins
1441 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1442 LOGICAL :: do_ppl, gapw, gapw_composite_direct_ao, &
1443 gapw_composite_reference, gapw_xc, &
1444 lrigpw, rigpw, use_work_v_rspace
1445 REAL(kind=dp) :: csign, dvol, fadm
1446 TYPE(admm_type), POINTER :: admm_env
1447 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1448 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao, rho_ao_nokp, smat
1449 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, &
1450 matrix_ks_aux_fit_dft, rho_ao_aux, &
1451 rho_ao_kp
1452 TYPE(dft_control_type), POINTER :: dft_control
1453 TYPE(kpoint_type), POINTER :: kpoints
1454 TYPE(lri_density_type), POINTER :: lri_density
1455 TYPE(lri_environment_type), POINTER :: lri_env
1456 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
1457 TYPE(mp_para_env_type), POINTER :: para_env
1458 TYPE(pw_env_type), POINTER :: pw_env
1459 TYPE(pw_poisson_type), POINTER :: poisson_env
1460 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1461 TYPE(pw_r3d_rs_type), POINTER :: v_rspace, v_rspace_used, vee
1462 TYPE(pw_r3d_rs_type), TARGET :: v_rspace_work
1463 TYPE(qs_ks_env_type), POINTER :: ks_env
1464 TYPE(qs_rho_type), POINTER :: rho_aux_fit
1465 TYPE(section_vals_type), POINTER :: input, xc_section
1466 TYPE(task_list_type), POINTER :: task_list
1467
1468 CALL timeset(routinen, handle)
1469
1470 NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, &
1471 v_rspace, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, &
1472 ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, &
1473 rho_ao_nokp, ks_env, admm_env, task_list, v_rspace_used, input, xc_section)
1474
1475 CALL get_qs_env(qs_env, &
1476 dft_control=dft_control, &
1477 input=input, &
1478 pw_env=pw_env, &
1479 v_hartree_rspace=v_rspace, &
1480 vee=vee)
1481
1482 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1483 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
1484 gapw = dft_control%qs_control%gapw
1485 gapw_xc = dft_control%qs_control%gapw_xc
1486 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1487 gapw_composite_reference = native_skala_gapw_composite_reference(xc_section) .AND. &
1488 (gapw .OR. gapw_xc)
1489 gapw_composite_direct_ao = gapw_composite_reference .AND. &
1490 native_skala_gapw_composite_direct_ao(xc_section)
1491 do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
1492
1493 rigpw = dft_control%qs_control%rigpw
1494 lrigpw = dft_control%qs_control%lrigpw
1495 IF (lrigpw .OR. rigpw) THEN
1496 CALL get_qs_env(qs_env, &
1497 lri_env=lri_env, &
1498 lri_density=lri_density, &
1499 atomic_kind_set=atomic_kind_set)
1500 END IF
1501
1502 nspins = dft_control%nspins
1503
1504 ! sum up potentials and integrate
1505 IF (ASSOCIATED(v_rspace_new)) THEN
1506 DO ispin = 1, nspins
1507 IF (gapw_composite_reference) THEN
1508 ! The direct AO diagnostic uses the full ORB basis. The reconstructed path uses
1509 ! the soft basis here and adds the hard-minus-soft adjoint through update_ks_atom.
1510 CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
1511 rho_ao => rho_ao_kp(ispin, :)
1512 ksmat => ks_matrix(ispin, :)
1513 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1514 pmat_kp=rho_ao, hmat_kp=ksmat, &
1515 qs_env=qs_env, &
1516 calculate_forces=calculate_forces, &
1517 gapw=gapw .AND. .NOT. gapw_composite_direct_ao)
1518 CALL pw_copy(v_rspace, v_rspace_new(ispin))
1519 ELSE IF (gapw_xc) THEN
1520 ! SIC not implemented (or at least not tested)
1521 cpassert(dft_control%sic_method_id == sic_none)
1522 !Only the xc potential, because it has to be integrated with the soft basis
1523 CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
1524
1525 ! add the xc part due to v_rspace soft
1526 rho_ao => rho_ao_kp(ispin, :)
1527 ksmat => ks_matrix(ispin, :)
1528 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1529 pmat_kp=rho_ao, hmat_kp=ksmat, &
1530 qs_env=qs_env, &
1531 calculate_forces=calculate_forces, &
1532 gapw=gapw_xc)
1533
1534 ! Now the Hartree potential to be integrated with the full basis
1535 CALL pw_copy(v_rspace, v_rspace_new(ispin))
1536 ELSE
1537 ! Add v_hartree + v_xc = v_rspace_new
1538 CALL pw_axpy(v_rspace, v_rspace_new(ispin), 1.0_dp, v_rspace_new(ispin)%pw_grid%dvol)
1539 END IF ! gapw_xc
1540 IF (dft_control%qs_control%ddapc_explicit_potential) THEN
1541 IF (dft_control%qs_control%ddapc_restraint_is_spin) THEN
1542 IF (ispin == 1) THEN
1543 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1544 ELSE
1545 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), -1.0_dp)
1546 END IF
1547 ELSE
1548 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1549 END IF
1550 END IF
1551 ! CDFT constraint contribution
1552 IF (dft_control%qs_control%cdft) THEN
1553 DO igroup = 1, SIZE(cdft_control%group)
1554 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1555 CASE (cdft_charge_constraint)
1556 csign = 1.0_dp
1557 CASE (cdft_magnetization_constraint)
1558 IF (ispin == 1) THEN
1559 csign = 1.0_dp
1560 ELSE
1561 csign = -1.0_dp
1562 END IF
1563 CASE (cdft_alpha_constraint)
1564 csign = 1.0_dp
1565 IF (ispin == 2) cycle
1566 CASE (cdft_beta_constraint)
1567 csign = 1.0_dp
1568 IF (ispin == 1) cycle
1569 CASE DEFAULT
1570 cpabort("Unknown constraint type.")
1571 END SELECT
1572 CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_new(ispin), &
1573 csign*cdft_control%strength(igroup))
1574 END DO
1575 END IF
1576 ! functional derivative of the Hartree energy wrt the density in the presence of dielectric
1577 ! (vhartree + v_eps); v_eps is nonzero only if the dielectric constant is defind as a function
1578 ! of the charge density
1579 IF (poisson_env%parameters%solver == pw_poisson_implicit) THEN
1580 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1581 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_new(ispin), dvol)
1582 END IF
1583 ! Add SCCS contribution
1584 IF (dft_control%do_sccs) THEN
1585 CALL pw_axpy(v_sccs_rspace, v_rspace_new(ispin))
1586 END IF
1587 ! External electrostatic potential
1588 IF (dft_control%apply_external_potential) THEN
1589 CALL qmmm_modify_hartree_pot(v_hartree=v_rspace_new(ispin), &
1590 v_qmmm=vee, scale=-1.0_dp)
1591 END IF
1592 IF (do_ppl) THEN
1593 cpassert(.NOT. gapw)
1594 CALL pw_axpy(vppl_rspace, v_rspace_new(ispin), vppl_rspace%pw_grid%dvol)
1595 END IF
1596 ! the electrostatic sic contribution
1597 SELECT CASE (dft_control%sic_method_id)
1598 CASE (sic_none)
1599 !
1600 CASE (sic_mauri_us, sic_mauri_spz)
1601 IF (ispin == 1) THEN
1602 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1603 ELSE
1604 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), 1.0_dp)
1605 END IF
1606 CASE (sic_ad)
1607 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1608 CASE (sic_eo)
1609 ! NOTHING TO BE DONE
1610 END SELECT
1611 ! DFT embedding
1612 IF (dft_control%apply_embed_pot) THEN
1613 CALL pw_axpy(v_rspace_embed(ispin), v_rspace_new(ispin), v_rspace_embed(ispin)%pw_grid%dvol)
1614 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1615 END IF
1616 IF (lrigpw) THEN
1617 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1618 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1619 DO ikind = 1, nkind
1620 lri_v_int(ikind)%v_int = 0.0_dp
1621 END DO
1622 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1623 lri_v_int, calculate_forces, "LRI_AUX")
1624 DO ikind = 1, nkind
1625 CALL para_env%sum(lri_v_int(ikind)%v_int)
1626 END DO
1627 IF (lri_env%exact_1c_terms) THEN
1628 rho_ao => my_rho(ispin, :)
1629 ksmat => ks_matrix(ispin, :)
1630 CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), ksmat(1)%matrix, &
1631 rho_ao(1)%matrix, qs_env, &
1632 calculate_forces, "ORB")
1633 END IF
1634 IF (lri_env%ppl_ri) THEN
1635 CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1636 END IF
1637 ELSE IF (rigpw) THEN
1638 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1639 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1640 DO ikind = 1, nkind
1641 lri_v_int(ikind)%v_int = 0.0_dp
1642 END DO
1643 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1644 lri_v_int, calculate_forces, "RI_HXC")
1645 DO ikind = 1, nkind
1646 CALL para_env%sum(lri_v_int(ikind)%v_int)
1647 END DO
1648 ELSE
1649 rho_ao => my_rho(ispin, :)
1650 ksmat => ks_matrix(ispin, :)
1651 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1652 pmat_kp=rho_ao, hmat_kp=ksmat, &
1653 qs_env=qs_env, &
1654 calculate_forces=calculate_forces, &
1655 gapw=gapw .AND. .NOT. gapw_composite_direct_ao)
1656 END IF
1657 CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
1658 END DO ! ispin
1659
1660 SELECT CASE (dft_control%sic_method_id)
1661 CASE (sic_none)
1662 CASE (sic_mauri_us, sic_mauri_spz, sic_ad)
1663 CALL auxbas_pw_pool%give_back_pw(v_sic_rspace)
1664 DEALLOCATE (v_sic_rspace)
1665 END SELECT
1666 DEALLOCATE (v_rspace_new)
1667
1668 ELSE
1669 ! not implemented (or at least not tested)
1670 cpassert(dft_control%sic_method_id == sic_none)
1671 cpassert(.NOT. dft_control%qs_control%ddapc_restraint_is_spin)
1672 DO ispin = 1, nspins
1673 use_work_v_rspace = dft_control%qs_control%cdft
1674 IF (use_work_v_rspace) THEN
1675 CALL auxbas_pw_pool%create_pw(v_rspace_work)
1676 CALL pw_copy(v_rspace, v_rspace_work)
1677 v_rspace_used => v_rspace_work
1678 ELSE
1679 v_rspace_used => v_rspace
1680 END IF
1681 ! CDFT constraint contribution
1682 IF (dft_control%qs_control%cdft) THEN
1683 DO igroup = 1, SIZE(cdft_control%group)
1684 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1685 CASE (cdft_charge_constraint)
1686 csign = 1.0_dp
1687 CASE (cdft_magnetization_constraint)
1688 IF (ispin == 1) THEN
1689 csign = 1.0_dp
1690 ELSE
1691 csign = -1.0_dp
1692 END IF
1693 CASE (cdft_alpha_constraint)
1694 csign = 1.0_dp
1695 IF (ispin == 2) cycle
1696 CASE (cdft_beta_constraint)
1697 csign = 1.0_dp
1698 IF (ispin == 1) cycle
1699 CASE DEFAULT
1700 cpabort("Unknown constraint type.")
1701 END SELECT
1702 CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_used, &
1703 csign*cdft_control%strength(igroup))
1704 END DO
1705 END IF
1706 ! extra contribution attributed to the dependency of the dielectric constant to the charge density
1707 IF (poisson_env%parameters%solver == pw_poisson_implicit) THEN
1708 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1709 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_used, dvol)
1710 END IF
1711 ! Add SCCS contribution
1712 IF (dft_control%do_sccs) THEN
1713 CALL pw_axpy(v_sccs_rspace, v_rspace_used)
1714 END IF
1715 ! DFT embedding
1716 IF (dft_control%apply_embed_pot) THEN
1717 CALL pw_axpy(v_rspace_embed(ispin), v_rspace_used, v_rspace_embed(ispin)%pw_grid%dvol)
1718 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1719 END IF
1720 IF (lrigpw) THEN
1721 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1722 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1723 DO ikind = 1, nkind
1724 lri_v_int(ikind)%v_int = 0.0_dp
1725 END DO
1726 CALL integrate_v_rspace_one_center(v_rspace_used, qs_env, &
1727 lri_v_int, calculate_forces, "LRI_AUX")
1728 DO ikind = 1, nkind
1729 CALL para_env%sum(lri_v_int(ikind)%v_int)
1730 END DO
1731 IF (lri_env%exact_1c_terms) THEN
1732 rho_ao => my_rho(ispin, :)
1733 ksmat => ks_matrix(ispin, :)
1734 CALL integrate_v_rspace_diagonal(v_rspace_used, ksmat(1)%matrix, &
1735 rho_ao(1)%matrix, qs_env, &
1736 calculate_forces, "ORB")
1737 END IF
1738 IF (lri_env%ppl_ri) THEN
1739 CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1740 END IF
1741 ELSE IF (rigpw) THEN
1742 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1743 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1744 DO ikind = 1, nkind
1745 lri_v_int(ikind)%v_int = 0.0_dp
1746 END DO
1747 CALL integrate_v_rspace_one_center(v_rspace_used, qs_env, &
1748 lri_v_int, calculate_forces, "RI_HXC")
1749 DO ikind = 1, nkind
1750 CALL para_env%sum(lri_v_int(ikind)%v_int)
1751 END DO
1752 ELSE
1753 rho_ao => my_rho(ispin, :)
1754 ksmat => ks_matrix(ispin, :)
1755 CALL integrate_v_rspace(v_rspace=v_rspace_used, &
1756 pmat_kp=rho_ao, &
1757 hmat_kp=ksmat, &
1758 qs_env=qs_env, &
1759 calculate_forces=calculate_forces, &
1760 gapw=gapw)
1761 END IF
1762 IF (use_work_v_rspace) CALL auxbas_pw_pool%give_back_pw(v_rspace_work)
1763 END DO
1764 END IF ! ASSOCIATED(v_rspace_new)
1765
1766 ! **** LRIGPW: KS matrix from integrated potential
1767 IF (lrigpw) THEN
1768 CALL get_qs_env(qs_env, ks_env=ks_env)
1769 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1770 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1771 DO ispin = 1, nspins
1772 ksmat => ks_matrix(ispin, :)
1773 CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set, &
1774 cell_to_index=cell_to_index)
1775 END DO
1776 IF (calculate_forces) THEN
1777 CALL calculate_lri_forces(lri_env, lri_density, qs_env, rho_ao_kp, atomic_kind_set)
1778 END IF
1779 ELSE IF (rigpw) THEN
1780 CALL get_qs_env(qs_env, matrix_s=smat)
1781 DO ispin = 1, nspins
1782 CALL calculate_ri_ks_matrix(lri_env, lri_v_int, ks_matrix(ispin, 1)%matrix, &
1783 smat(1)%matrix, atomic_kind_set, ispin)
1784 END DO
1785 IF (calculate_forces) THEN
1786 rho_ao_nokp => rho_ao_kp(:, 1)
1787 CALL calculate_ri_forces(lri_env, lri_density, qs_env, rho_ao_nokp, atomic_kind_set)
1788 END IF
1789 END IF
1790
1791 IF (ASSOCIATED(v_tau_rspace)) THEN
1792 IF (lrigpw .OR. rigpw) THEN
1793 cpabort("LRIGPW/RIGPW not implemented for meta-GGAs")
1794 END IF
1795 DO ispin = 1, nspins
1796 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1797
1798 rho_ao => rho_ao_kp(ispin, :)
1799 ksmat => ks_matrix(ispin, :)
1800 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1801 pmat_kp=rho_ao, hmat_kp=ksmat, &
1802 qs_env=qs_env, &
1803 calculate_forces=calculate_forces, compute_tau=.true., &
1804 gapw=(gapw .OR. gapw_xc) .AND. &
1805 .NOT. gapw_composite_direct_ao)
1806 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1807 END DO
1808 DEALLOCATE (v_tau_rspace)
1809 END IF
1810
1811 ! Add contributions from ADMM if requested
1812 IF (dft_control%do_admm) THEN
1813 CALL get_qs_env(qs_env, admm_env=admm_env)
1814 CALL get_admm_env(admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
1815 matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft)
1816 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
1817 IF (ASSOCIATED(v_rspace_new_aux_fit)) THEN
1818 DO ispin = 1, nspins
1819 ! Calculate the xc potential
1820 CALL pw_scale(v_rspace_new_aux_fit(ispin), v_rspace_new_aux_fit(ispin)%pw_grid%dvol)
1821
1822 ! set matrix_ks_aux_fit_dft = matrix_ks_aux_fit(k_HF)
1823 DO img = 1, dft_control%nimages
1824 CALL dbcsr_copy(matrix_ks_aux_fit_dft(ispin, img)%matrix, matrix_ks_aux_fit(ispin, img)%matrix, &
1825 name="DFT exch. part of matrix_ks_aux_fit")
1826 END DO
1827
1828 ! Add potential to ks_matrix aux_fit, skip integration if no DFT correction
1829
1830 IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
1831
1832 !GPW by default. IF GAPW, then take relevant task list and basis
1833 CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
1834 basis_type = "AUX_FIT"
1835 IF (admm_env%do_gapw) THEN
1836 task_list => admm_env%admm_gapw_env%task_list
1837 basis_type = "AUX_FIT_SOFT"
1838 END IF
1839 fadm = 1.0_dp
1840 ! Calculate bare scaling of force according to Merlot, 1. IF: ADMMP, 2. IF: ADMMS,
1841 IF (admm_env%do_admmp) THEN
1842 fadm = admm_env%gsi(ispin)**2
1843 ELSE IF (admm_env%do_admms) THEN
1844 fadm = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)
1845 END IF
1846
1847 rho_ao => rho_ao_aux(ispin, :)
1848 ksmat => matrix_ks_aux_fit(ispin, :)
1849
1850 CALL integrate_v_rspace(v_rspace=v_rspace_new_aux_fit(ispin), &
1851 pmat_kp=rho_ao, &
1852 hmat_kp=ksmat, &
1853 qs_env=qs_env, &
1854 calculate_forces=calculate_forces, &
1855 force_adm=fadm, &
1856 gapw=.false., & !even if actual GAPW calculation, want to use AUX_FIT_SOFT
1857 basis_type=basis_type, &
1858 task_list_external=task_list)
1859 END IF
1860
1861 ! matrix_ks_aux_fit_dft(x_DFT)=matrix_ks_aux_fit_dft(old,k_HF)-matrix_ks_aux_fit(k_HF-x_DFT)
1862 DO img = 1, dft_control%nimages
1863 CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, &
1864 matrix_ks_aux_fit(ispin, img)%matrix, 1.0_dp, -1.0_dp)
1865 END DO
1866
1867 CALL auxbas_pw_pool%give_back_pw(v_rspace_new_aux_fit(ispin))
1868 END DO
1869 DEALLOCATE (v_rspace_new_aux_fit)
1870 END IF
1871 ! Clean up v_tau_rspace_aux_fit, which is actually not needed
1872 IF (ASSOCIATED(v_tau_rspace_aux_fit)) THEN
1873 DO ispin = 1, nspins
1874 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace_aux_fit(ispin))
1875 END DO
1876 DEALLOCATE (v_tau_rspace_aux_fit)
1877 END IF
1878 END IF
1879
1880 IF (dft_control%apply_embed_pot) DEALLOCATE (v_rspace_embed)
1881
1882 CALL timestop(handle)
1883
1884 END SUBROUTINE sum_up_and_integrate
1885
1886!**************************************************************************
1887!> \brief Calculate the ZMP potential and energy as in Zhao, Morrison Parr
1888!> PRA 50i, 2138 (1994)
1889!> V_c^\lambda defined as int_rho-rho_0/r-r' or rho-rho_0 times a Lagrange
1890!> multiplier, plus Fermi-Amaldi potential that should give the V_xc in the
1891!> limit \lambda --> \infty
1892!>
1893!> \param qs_env ...
1894!> \param v_rspace_new ...
1895!> \param rho ...
1896!> \param exc ...
1897!> \author D. Varsano [daniele.varsano@nano.cnr.it]
1898! **************************************************************************************************
1899 SUBROUTINE calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
1900
1901 TYPE(qs_environment_type), POINTER :: qs_env
1902 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new
1903 TYPE(qs_rho_type), POINTER :: rho
1904 REAL(kind=dp) :: exc
1905
1906 CHARACTER(*), PARAMETER :: routinen = 'calculate_zmp_potential'
1907
1908 INTEGER :: handle, my_val, nelectron, nspins
1909 INTEGER, DIMENSION(2) :: nelectron_spin
1910 LOGICAL :: do_zmp_read, fermi_amaldi
1911 REAL(kind=dp) :: lambda
1912 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_ext_r
1913 TYPE(dft_control_type), POINTER :: dft_control
1914 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ext_g, rho_g
1915 TYPE(pw_env_type), POINTER :: pw_env
1916 TYPE(pw_poisson_type), POINTER :: poisson_env
1917 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1918 TYPE(pw_r3d_rs_type) :: v_xc_rspace
1919 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1920 TYPE(qs_ks_env_type), POINTER :: ks_env
1921 TYPE(section_vals_type), POINTER :: ext_den_section, input
1922
1923!, v_h_gspace, &
1924
1925 CALL timeset(routinen, handle)
1926 NULLIFY (auxbas_pw_pool)
1927 NULLIFY (pw_env)
1928 NULLIFY (poisson_env)
1929 NULLIFY (v_rspace_new)
1930 NULLIFY (dft_control)
1931 NULLIFY (rho_r, rho_g, tot_rho_ext_r, rho_ext_g)
1932 CALL get_qs_env(qs_env=qs_env, &
1933 pw_env=pw_env, &
1934 ks_env=ks_env, &
1935 rho=rho, &
1936 input=input, &
1937 nelectron_spin=nelectron_spin, &
1938 dft_control=dft_control)
1939 CALL pw_env_get(pw_env=pw_env, &
1940 auxbas_pw_pool=auxbas_pw_pool, &
1941 poisson_env=poisson_env)
1942 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
1943 nspins = 1
1944 ALLOCATE (v_rspace_new(nspins))
1945 CALL auxbas_pw_pool%create_pw(pw=v_rspace_new(1))
1946 CALL auxbas_pw_pool%create_pw(pw=v_xc_rspace)
1947
1948 CALL pw_zero(v_rspace_new(1))
1949 do_zmp_read = dft_control%apply_external_vxc
1950 IF (do_zmp_read) THEN
1951 CALL pw_copy(qs_env%external_vxc, v_rspace_new(1))
1952 exc = accurate_dot_product(v_rspace_new(1)%array, rho_r(1)%array)* &
1953 v_rspace_new(1)%pw_grid%dvol
1954 ELSE
1955 block
1956 REAL(kind=dp) :: factor
1957 TYPE(pw_c1d_gs_type) :: rho_eff_gspace, v_xc_gspace
1958 CALL auxbas_pw_pool%create_pw(pw=rho_eff_gspace)
1959 CALL auxbas_pw_pool%create_pw(pw=v_xc_gspace)
1960 CALL pw_zero(rho_eff_gspace)
1961 CALL pw_zero(v_xc_gspace)
1962 CALL pw_zero(v_xc_rspace)
1963 factor = pw_integrate_function(rho_g(1))
1964 CALL qs_rho_get(qs_env%rho_external, &
1965 rho_g=rho_ext_g, &
1966 tot_rho_r=tot_rho_ext_r)
1967 factor = tot_rho_ext_r(1)/factor
1968
1969 CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1970 CALL pw_axpy(rho_ext_g(1), rho_eff_gspace, alpha=-1.0_dp)
1971 ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
1972 CALL section_vals_val_get(ext_den_section, "LAMBDA", r_val=lambda)
1973 CALL section_vals_val_get(ext_den_section, "ZMP_CONSTRAINT", i_val=my_val)
1974 CALL section_vals_val_get(ext_den_section, "FERMI_AMALDI", l_val=fermi_amaldi)
1975
1976 CALL pw_scale(rho_eff_gspace, a=lambda)
1977 nelectron = nelectron_spin(1)
1978 factor = -1.0_dp/nelectron
1979 CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1980
1981 CALL pw_poisson_solve(poisson_env, rho_eff_gspace, vhartree=v_xc_gspace)
1982 CALL pw_transfer(v_xc_gspace, v_rspace_new(1))
1983 CALL pw_copy(v_rspace_new(1), v_xc_rspace)
1984
1985 exc = 0.0_dp
1986 exc = pw_integral_ab(v_rspace_new(1), rho_r(1))
1987
1988!Note that this is not the xc energy but \int(\rho*v_xc)
1989!Vxc---> v_rspace_new
1990!Exc---> energy%exc
1991 CALL auxbas_pw_pool%give_back_pw(rho_eff_gspace)
1992 CALL auxbas_pw_pool%give_back_pw(v_xc_gspace)
1993 END block
1994 END IF
1995
1996 CALL auxbas_pw_pool%give_back_pw(v_xc_rspace)
1997
1998 CALL timestop(handle)
1999
2000 END SUBROUTINE calculate_zmp_potential
2001
2002! **************************************************************************************************
2003!> \brief ...
2004!> \param qs_env ...
2005!> \param rho ...
2006!> \param v_rspace_embed ...
2007!> \param dft_control ...
2008!> \param embed_corr ...
2009!> \param just_energy ...
2010! **************************************************************************************************
2011 SUBROUTINE get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, &
2012 just_energy)
2013 TYPE(qs_environment_type), POINTER :: qs_env
2014 TYPE(qs_rho_type), POINTER :: rho
2015 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_embed
2016 TYPE(dft_control_type), POINTER :: dft_control
2017 REAL(kind=dp) :: embed_corr
2018 LOGICAL :: just_energy
2019
2020 CHARACTER(*), PARAMETER :: routinen = 'get_embed_potential_energy'
2021
2022 INTEGER :: handle, ispin
2023 REAL(kind=dp) :: embed_corr_local
2024 TYPE(pw_env_type), POINTER :: pw_env
2025 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2026 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
2027
2028 CALL timeset(routinen, handle)
2029
2030 NULLIFY (auxbas_pw_pool)
2031 NULLIFY (pw_env)
2032 NULLIFY (rho_r)
2033 CALL get_qs_env(qs_env=qs_env, &
2034 pw_env=pw_env, &
2035 rho=rho)
2036 CALL pw_env_get(pw_env=pw_env, &
2037 auxbas_pw_pool=auxbas_pw_pool)
2038 CALL qs_rho_get(rho, rho_r=rho_r)
2039 ALLOCATE (v_rspace_embed(dft_control%nspins))
2040
2041 embed_corr = 0.0_dp
2042
2043 DO ispin = 1, dft_control%nspins
2044 CALL auxbas_pw_pool%create_pw(pw=v_rspace_embed(ispin))
2045 CALL pw_zero(v_rspace_embed(ispin))
2046
2047 CALL pw_copy(qs_env%embed_pot, v_rspace_embed(ispin))
2048 embed_corr_local = 0.0_dp
2049
2050 ! Spin embedding potential in open-shell case
2051 IF (dft_control%nspins == 2) THEN
2052 IF (ispin == 1) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), 1.0_dp)
2053 IF (ispin == 2) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), -1.0_dp)
2054 END IF
2055 ! Integrate the density*potential
2056 embed_corr_local = pw_integral_ab(v_rspace_embed(ispin), rho_r(ispin))
2057
2058 embed_corr = embed_corr + embed_corr_local
2059
2060 END DO
2061
2062 ! If only energy requiested we delete the potential
2063 IF (just_energy) THEN
2064 DO ispin = 1, dft_control%nspins
2065 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
2066 END DO
2067 DEALLOCATE (v_rspace_embed)
2068 END IF
2069
2070 CALL timestop(handle)
2071
2072 END SUBROUTINE get_embed_potential_energy
2073
2074END MODULE qs_ks_utils
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:95
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:103
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
cp_fm_types::cp_fm_to_fm
Definition cp_fm_types.F:83
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
kahan_sum::accurate_dot_product
Definition kahan_sum.F:52
kahan_sum::accurate_sum
Definition kahan_sum.F:41
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_integrate_function
Definition pw_methods.F:116
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:599
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
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_release_p
subroutine, public dbcsr_release_p(matrix)
...
Definition cp_dbcsr_api.F:222
cp_dbcsr_api::dbcsr_scale
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
Definition cp_dbcsr_api.F:1188
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_contrib
Definition cp_dbcsr_contrib.F:8
cp_dbcsr_contrib::dbcsr_dot
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
Definition cp_dbcsr_contrib.F:367
cp_dbcsr_contrib::dbcsr_scale_by_vector
subroutine, public dbcsr_scale_by_vector(matrix, alpha, side)
Scales the rows/columns of given matrix.
Definition cp_dbcsr_contrib.F:486
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition cp_dbcsr_operations.F:799
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm, plan)
Copy a DBCSR matrix to a BLACS matrix.
Definition cp_dbcsr_operations.F:232
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:944
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition cp_dbcsr_operations.F:129
cp_ddapc
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition cp_ddapc.F:15
cp_ddapc::cp_ddapc_apply_cd
subroutine, public cp_ddapc_apply_cd(qs_env, rho_tot_gspace, energy, v_hartree_gspace, calculate_forces, itype_of_density)
Routine to couple/decouple periodic images with the Bloechl scheme.
Definition cp_ddapc.F:226
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:132
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:355
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
hfx_admm_utils
Utilities for hfx and admm methods.
Definition hfx_admm_utils.F:17
hfx_admm_utils::tddft_hfx_matrix
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data, external_para_env)
Add the hfx contributions to the Hamiltonian.
Definition hfx_admm_utils.F:2489
hfx_derivatives
Routines to calculate derivatives with respect to basis function origin.
Definition hfx_derivatives.F:14
hfx_derivatives::derivatives_four_center
subroutine, public derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, irep, use_virial, adiabatic_rescale_factor, resp_only, external_x_data, nspins)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Definition hfx_derivatives.F:117
hfx_types
Types and set/get functions for HFX.
Definition hfx_types.F:16
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::sic_list_unpaired
integer, parameter, public sic_list_unpaired
Definition input_constants.F:615
input_constants::sic_mauri_spz
integer, parameter, public sic_mauri_spz
Definition input_constants.F:610
input_constants::cdft_beta_constraint
integer, parameter, public cdft_beta_constraint
Definition input_constants.F:820
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition input_constants.F:820
input_constants::sic_list_all
integer, parameter, public sic_list_all
Definition input_constants.F:615
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition input_constants.F:820
input_constants::sic_eo
integer, parameter, public sic_eo
Definition input_constants.F:610
input_constants::do_ppl_grid
integer, parameter, public do_ppl_grid
Definition input_constants.F:294
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition input_constants.F:921
input_constants::sic_mauri_us
integer, parameter, public sic_mauri_us
Definition input_constants.F:610
input_constants::sic_none
integer, parameter, public sic_none
Definition input_constants.F:610
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition input_constants.F:820
input_constants::sic_ad
integer, parameter, public sic_ad
Definition input_constants.F:610
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
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
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:587
lri_environment_methods
Calculates integral matrices for LRIGPW method lri : local resolution of the identity.
Definition lri_environment_methods.F:18
lri_environment_methods::v_int_ppl_update
subroutine, public v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
...
Definition lri_environment_methods.F:1189
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition lri_environment_types.F:18
lri_forces
Calculates forces for LRIGPW method lri : local resolution of the identity.
Definition lri_forces.F:16
lri_forces::calculate_lri_forces
subroutine, public calculate_lri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the lri forces
Definition lri_forces.F:81
lri_forces::calculate_ri_forces
subroutine, public calculate_ri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the ri forces
Definition lri_forces.F:673
lri_ks_methods
routines that build the Kohn-Sham matrix for the LRIGPW and xc parts
Definition lri_ks_methods.F:15
lri_ks_methods::calculate_lri_ks_matrix
subroutine, public calculate_lri_ks_matrix(lri_env, lri_v_int, h_matrix, atomic_kind_set, cell_to_index)
update of LRIGPW KS matrix
Definition lri_ks_methods.F:65
lri_ks_methods::calculate_ri_ks_matrix
subroutine, public calculate_ri_ks_matrix(lri_env, lri_v_int, h_matrix, s_matrix, atomic_kind_set, ispin)
update of RIGPW KS matrix
Definition lri_ks_methods.F:285
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::abnormal_value
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Definition mathlib.F:159
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
ps_implicit_types
Types containing essential information for running implicit (iterative) Poisson solver.
Definition ps_implicit_types.F:15
ps_implicit_types::neumann_bc
integer, parameter, public neumann_bc
Definition ps_implicit_types.F:32
ps_implicit_types::mixed_bc
integer, parameter, public mixed_bc
Definition ps_implicit_types.F:32
ps_implicit_types::mixed_periodic_bc
integer, parameter, public mixed_periodic_bc
Definition ps_implicit_types.F:32
ps_implicit_types::periodic_bc
integer, parameter, public periodic_bc
Definition ps_implicit_types.F:32
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_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_poisson_types::pw_poisson_implicit
integer, parameter, public pw_poisson_implicit
Definition pw_poisson_types.F:67
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_cdft_types
Defines CDFT control structures.
Definition qs_cdft_types.F:14
qs_charges_types
container for information about total charges on the grids
Definition qs_charges_types.F:14
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition qs_collocate_density.F:1720
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:545
qs_force_types
Definition qs_force_types.F:13
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_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, do_mtlr_present, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
Definition qs_kind_types.F:984
qs_ks_qmmm_methods
Definition qs_ks_qmmm_methods.F:11
qs_ks_qmmm_methods::qmmm_modify_hartree_pot
subroutine, public qmmm_modify_hartree_pot(v_hartree, v_qmmm, scale)
Modify the hartree potential in order to include the QM/MM correction.
Definition qs_ks_qmmm_methods.F:205
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition qs_ks_types.F:343
qs_ks_utils
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_utils.F:22
qs_ks_utils::print_densities
subroutine, public print_densities(qs_env, rho)
...
Definition qs_ks_utils.F:958
qs_ks_utils::get_embed_potential_energy
subroutine, public get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, just_energy)
...
Definition qs_ks_utils.F:2013
qs_ks_utils::compute_matrix_vxc_kp
subroutine, public compute_matrix_vxc_kp(qs_env, v_rspace, matrix_vxc_kp, gapw_full_basis)
Build the XC potential matrix for k-point/image-resolved KS matrices.
Definition qs_ks_utils.F:1339
qs_ks_utils::low_spin_roks
subroutine, public low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, calculate_forces, auxbas_pw_pool)
do ROKS calculations yielding low spin states
Definition qs_ks_utils.F:153
qs_ks_utils::check_sum_energies
subroutine, public check_sum_energies(energy)
Check each term of the energy and sum to total.
Definition qs_ks_utils.F:891
qs_ks_utils::compute_matrix_vxc
subroutine, public compute_matrix_vxc(qs_env, v_rspace, matrix_vxc, gapw_full_basis)
compute matrix_vxc, defined via the potential created by qs_vxc_create ignores things like tau functi...
Definition qs_ks_utils.F:1282
qs_ks_utils::sum_up_and_integrate
subroutine, public sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, cdft_control, calculate_forces)
Sum up all potentials defined on the grid and integrate.
Definition qs_ks_utils.F:1422
qs_ks_utils::print_detailed_energy
subroutine, public print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
Print detailed energies.
Definition qs_ks_utils.F:1106
qs_ks_utils::calculate_zmp_potential
subroutine, public calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
Calculate the ZMP potential and energy as in Zhao, Morrison Parr PRA 50i, 2138 (1994) V_c^\lambda def...
Definition qs_ks_utils.F:1900
qs_ks_utils::calc_v_sic_rspace
subroutine, public calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, just_energy, calculate_forces, auxbas_pw_pool)
do sic calculations on the spin density
Definition qs_ks_utils.F:772
qs_ks_utils::sic_explicit_orbitals
subroutine, public sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, calculate_forces, auxbas_pw_pool)
do sic calculations on explicit orbitals
Definition qs_ks_utils.F:476
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, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:445
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
skala_gpw_functional
Experimental CP2K-native GPW real-space-grid path for SKALA TorchScript models.
Definition skala_gpw_functional.F:11
skala_gpw_functional::native_skala_gapw_composite_direct_ao
logical function, public native_skala_gapw_composite_direct_ao(xc_section)
Return true if the GAPW composite reference uses direct full-ORB collocation.
Definition skala_gpw_functional.F:258
skala_gpw_functional::native_skala_gapw_composite_reference
logical function, public native_skala_gapw_composite_reference(xc_section)
Return true if native SKALA should use the full GAPW ORB density on one common grid.
Definition skala_gpw_functional.F:179
task_list_types
types for task lists
Definition task_list_types.F:14
virial_types
Definition virial_types.F:13
xc
Exchange and Correlation functional calculations.
Definition xc.F:17
xc::xc_exc_calc
real(kind=dp) function, public xc_exc_calc(rho_r, rho_g, tau, xc_section, weights, pw_pool)
calculates just the exchange and correlation energy (no vxc)
Definition xc.F:800
xc::xc_vxc_pw_create
subroutine, public xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, weights, pw_pool, compute_virial, virial_xc, exc_r)
Exchange and Correlation functional calculations.
Definition xc.F:483
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:60
cp_control_types::dft_control_type
Definition cp_control_types.F:780
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
hfx_types::hfx_type
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:514
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:198
lri_environment_types::lri_density_type
Definition lri_environment_types.F:359
lri_environment_types::lri_environment_type
Definition lri_environment_types.F:271
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:332
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_cdft_types::cdft_control_type
Definition qs_cdft_types.F:219
qs_charges_types::qs_charges_type
Container for information about total charges on the grids.
Definition qs_charges_types.F:41
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
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:190
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
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