Line data Source code
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 :
22 : MODULE qs_ks_utils
23 : USE admm_types, ONLY: admm_type,&
24 : get_admm_env
25 : USE atomic_kind_types, ONLY: atomic_kind_type,&
26 : get_atomic_kind
27 : USE cell_types, ONLY: cell_type
28 : USE cp_control_types, ONLY: dft_control_type
29 : USE cp_dbcsr_api, ONLY: &
30 : dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_info, dbcsr_init_p, &
31 : dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_scale, dbcsr_set, dbcsr_type
32 : USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
33 : dbcsr_scale_by_vector
34 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
35 : copy_fm_to_dbcsr,&
36 : cp_dbcsr_plus_fm_fm_t,&
37 : cp_dbcsr_sm_fm_multiply,&
38 : dbcsr_allocate_matrix_set,&
39 : dbcsr_deallocate_matrix_set
40 : USE cp_ddapc, ONLY: cp_ddapc_apply_CD
41 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
42 : cp_fm_struct_release,&
43 : cp_fm_struct_type
44 : USE cp_fm_types, ONLY: cp_fm_create,&
45 : cp_fm_get_info,&
46 : cp_fm_release,&
47 : cp_fm_set_all,&
48 : cp_fm_to_fm,&
49 : cp_fm_type
50 : USE cp_log_handling, ONLY: cp_get_default_logger,&
51 : cp_logger_type,&
52 : cp_to_string
53 : USE cp_output_handling, ONLY: cp_p_file,&
54 : cp_print_key_finished_output,&
55 : cp_print_key_should_output,&
56 : cp_print_key_unit_nr
57 : USE hfx_admm_utils, ONLY: tddft_hfx_matrix
58 : USE hfx_derivatives, ONLY: derivatives_four_center
59 : USE hfx_types, ONLY: hfx_type
60 : USE input_constants, ONLY: &
61 : cdft_alpha_constraint, cdft_beta_constraint, cdft_charge_constraint, &
62 : cdft_magnetization_constraint, do_admm_aux_exch_func_none, do_ppl_grid, sic_ad, sic_eo, &
63 : sic_list_all, sic_list_unpaired, sic_mauri_spz, sic_mauri_us, sic_none
64 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
65 : section_vals_type,&
66 : section_vals_val_get
67 : USE kahan_sum, ONLY: accurate_dot_product,&
68 : accurate_sum
69 : USE kinds, ONLY: default_string_length,&
70 : dp
71 : USE kpoint_types, ONLY: get_kpoint_info,&
72 : kpoint_type
73 : USE lri_environment_methods, ONLY: v_int_ppl_update
74 : USE lri_environment_types, ONLY: lri_density_type,&
75 : lri_environment_type,&
76 : lri_kind_type
77 : USE lri_forces, ONLY: calculate_lri_forces,&
78 : calculate_ri_forces
79 : USE lri_ks_methods, ONLY: calculate_lri_ks_matrix,&
80 : calculate_ri_ks_matrix
81 : USE mathlib, ONLY: abnormal_value
82 : USE message_passing, ONLY: mp_para_env_type
83 : USE ps_implicit_types, ONLY: MIXED_BC,&
84 : MIXED_PERIODIC_BC,&
85 : NEUMANN_BC,&
86 : PERIODIC_BC
87 : USE pw_env_types, ONLY: pw_env_get,&
88 : pw_env_type
89 : USE pw_methods, ONLY: pw_axpy,&
90 : pw_copy,&
91 : pw_integral_ab,&
92 : pw_integrate_function,&
93 : pw_scale,&
94 : pw_transfer,&
95 : pw_zero
96 : USE pw_poisson_methods, ONLY: pw_poisson_solve
97 : USE pw_poisson_types, ONLY: pw_poisson_implicit,&
98 : pw_poisson_type
99 : USE pw_pool_types, ONLY: pw_pool_type
100 : USE pw_types, ONLY: pw_c1d_gs_type,&
101 : pw_r3d_rs_type
102 : USE qs_cdft_types, ONLY: cdft_control_type
103 : USE qs_charges_types, ONLY: qs_charges_type
104 : USE qs_collocate_density, ONLY: calculate_rho_elec
105 : USE qs_energy_types, ONLY: qs_energy_type
106 : USE qs_environment_types, ONLY: get_qs_env,&
107 : qs_environment_type
108 : USE qs_force_types, ONLY: qs_force_type
109 : USE qs_integrate_potential, ONLY: integrate_v_rspace,&
110 : integrate_v_rspace_diagonal,&
111 : integrate_v_rspace_one_center
112 : USE qs_kind_types, ONLY: get_qs_kind,&
113 : get_qs_kind_set,&
114 : qs_kind_type
115 : USE qs_ks_atom, ONLY: update_ks_atom
116 : USE qs_ks_qmmm_methods, ONLY: qmmm_modify_hartree_pot
117 : USE qs_ks_types, ONLY: get_ks_env,&
118 : qs_ks_env_type
119 : USE qs_mo_types, ONLY: get_mo_set,&
120 : mo_set_type
121 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
122 : USE qs_oce_types, ONLY: oce_matrix_type
123 : USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
124 : calculate_rho_atom,&
125 : calculate_rho_atom_coeff
126 : USE qs_rho_atom_types, ONLY: deallocate_rho_atom_set,&
127 : rho_atom_type,&
128 : zero_rho_atom_integrals
129 : USE qs_rho_types, ONLY: qs_rho_get,&
130 : qs_rho_type
131 : USE qs_vxc_atom, ONLY: calculate_vxc_atom
132 : USE skala_gpw_functional, ONLY: native_skala_gapw_composite_direct_ao,&
133 : native_skala_gapw_composite_reference
134 : USE task_list_types, ONLY: task_list_type
135 : USE virial_types, ONLY: virial_type
136 : USE xc, ONLY: xc_exc_calc,&
137 : xc_vxc_pw_create
138 : USE xc_derivatives, ONLY: xc_functionals_get_needs
139 : USE xc_gauxc_functional, ONLY: apply_gauxc
140 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
141 : #include "./base/base_uses.f90"
142 :
143 : IMPLICIT NONE
144 :
145 : PRIVATE
146 :
147 : LOGICAL, PARAMETER :: debug_this_module = .TRUE.
148 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_utils'
149 :
150 : PUBLIC :: low_spin_roks, sic_explicit_orbitals, calc_v_sic_rspace, check_sum_energies, print_densities, &
151 : print_detailed_energy, compute_matrix_vxc, compute_matrix_vxc_kp, sum_up_and_integrate, &
152 : calculate_zmp_potential, get_embed_potential_energy
153 :
154 : CONTAINS
155 :
156 : ! **************************************************************************************************
157 : !> \brief do ROKS calculations yielding low spin states
158 : !> \param energy ...
159 : !> \param qs_env ...
160 : !> \param dft_control ...
161 : !> \param do_hfx ...
162 : !> \param just_energy ...
163 : !> \param calculate_forces ...
164 : !> \param auxbas_pw_pool ...
165 : ! **************************************************************************************************
166 136161 : SUBROUTINE low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
167 : calculate_forces, auxbas_pw_pool)
168 :
169 : TYPE(qs_energy_type), POINTER :: energy
170 : TYPE(qs_environment_type), POINTER :: qs_env
171 : TYPE(dft_control_type), POINTER :: dft_control
172 : LOGICAL, INTENT(IN) :: do_hfx, just_energy, calculate_forces
173 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
174 :
175 : CHARACTER(*), PARAMETER :: routineN = 'low_spin_roks'
176 :
177 : INTEGER :: handle, irep, ispin, iterm, k, k_alpha, &
178 : k_beta, n_rep, Nelectron, Nspin, Nterms
179 136161 : INTEGER, DIMENSION(:), POINTER :: ivec
180 136161 : INTEGER, DIMENSION(:, :, :), POINTER :: occupations
181 : LOGICAL :: compute_tau, compute_virial, gapw, &
182 : in_range, uniform_occupation, use_gauxc
183 : REAL(KIND=dp) :: ehfx, exc, exc1
184 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp
185 136161 : REAL(KIND=dp), DIMENSION(:), POINTER :: energy_scaling, rvec, scaling
186 : TYPE(cell_type), POINTER :: cell
187 136161 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_hfx, matrix_p, mdummy, &
188 136161 : mo_derivs, rho_ao
189 136161 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p2
190 : TYPE(dbcsr_type), POINTER :: dbcsr_deriv, fm_deriv, fm_scaled, &
191 : mo_coeff
192 136161 : TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
193 136161 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
194 : TYPE(mp_para_env_type), POINTER :: para_env
195 136161 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, tau_g
196 : TYPE(pw_env_type), POINTER :: pw_env
197 : TYPE(pw_pool_type), POINTER :: xc_pw_pool
198 : TYPE(pw_r3d_rs_type) :: work_v_rspace
199 136161 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc, vxc_tau
200 : TYPE(pw_r3d_rs_type), POINTER :: weights
201 : TYPE(qs_ks_env_type), POINTER :: ks_env
202 : TYPE(qs_rho_type), POINTER :: rho
203 136161 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
204 : TYPE(section_vals_type), POINTER :: hfx_section, input, &
205 : low_spin_roks_section, xc_section
206 : TYPE(virial_type), POINTER :: virial
207 : TYPE(xc_rho_cflags_type) :: needs
208 :
209 135663 : IF (.NOT. dft_control%low_spin_roks) RETURN
210 :
211 498 : CALL timeset(routineN, handle)
212 :
213 498 : NULLIFY (ks_env, rho_ao, rho_atom_set)
214 498 : gapw = dft_control%qs_control%gapw
215 498 : use_gauxc = gapw .AND. dft_control%use_gauxc
216 498 : compute_tau = .FALSE.
217 :
218 : ! Test for not compatible options
219 498 : IF (dft_control%qs_control%gapw_xc) THEN
220 0 : CALL cp_abort(__LOCATION__, "GAPW_XC is not compatible with low spin ROKS method.")
221 : END IF
222 498 : IF (dft_control%do_admm) THEN
223 0 : CALL cp_abort(__LOCATION__, "ADMM not compatible with low spin ROKS method.")
224 : END IF
225 498 : IF (dft_control%do_admm) THEN
226 0 : IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
227 : CALL cp_abort(__LOCATION__, "ADMM with XC correction functional "// &
228 0 : "not compatible with low spin ROKS method.")
229 : END IF
230 : END IF
231 498 : IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
232 : dft_control%qs_control%xtb) THEN
233 0 : CALL cp_abort(__LOCATION__, "SE/xTB/DFTB are not compatible with low spin ROKS method.")
234 : END IF
235 :
236 : CALL get_qs_env(qs_env, &
237 : ks_env=ks_env, &
238 : mo_derivs=mo_derivs, &
239 : mos=mo_array, &
240 : rho=rho, &
241 : pw_env=pw_env, &
242 : xcint_weights=weights, &
243 : input=input, &
244 : cell=cell, &
245 498 : virial=virial)
246 :
247 498 : CALL qs_rho_get(rho, rho_ao=rho_ao)
248 :
249 498 : compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
250 498 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
251 498 : hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
252 :
253 498 : IF (gapw) THEN
254 136 : IF (compute_virial) CPABORT("GAPW low spin ROKS stress is not implemented.")
255 136 : IF (calculate_forces) CPABORT("GAPW low spin ROKS forces are not implemented.")
256 136 : IF (.NOT. use_gauxc) THEN
257 : needs = xc_functionals_get_needs(section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL"), &
258 104 : lsd=.TRUE., calc_potential=.TRUE.)
259 104 : compute_tau = needs%tau .OR. needs%tau_spin
260 : END IF
261 : END IF
262 :
263 : ! The configuration corrections currently use unweighted grid integration.
264 498 : IF (ASSOCIATED(weights)) THEN
265 0 : CALL cp_abort(__LOCATION__, "No accurate xc integration possible.")
266 : END IF
267 : ! some assumptions need to be checked
268 : ! we have two spins
269 498 : CPASSERT(SIZE(mo_array, 1) == 2)
270 498 : Nspin = 2
271 : ! we want uniform occupations
272 498 : CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
273 498 : CPASSERT(uniform_occupation)
274 498 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation)
275 498 : CPASSERT(uniform_occupation)
276 498 : IF (do_hfx .AND. calculate_forces .AND. compute_virial) THEN
277 0 : CALL cp_abort(__LOCATION__, "ROKS virial with HFX not available.")
278 : END IF
279 :
280 498 : NULLIFY (dbcsr_deriv)
281 498 : CALL dbcsr_init_p(dbcsr_deriv)
282 498 : CALL dbcsr_copy(dbcsr_deriv, mo_derivs(1)%matrix)
283 498 : CALL dbcsr_set(dbcsr_deriv, 0.0_dp)
284 :
285 : ! basic info
286 498 : CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
287 498 : CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_alpha)
288 498 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff)
289 498 : CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_beta)
290 :
291 : ! read the input
292 498 : low_spin_roks_section => section_vals_get_subs_vals(input, "DFT%LOW_SPIN_ROKS")
293 :
294 498 : CALL section_vals_val_get(low_spin_roks_section, "ENERGY_SCALING", r_vals=rvec)
295 498 : Nterms = SIZE(rvec)
296 1494 : ALLOCATE (energy_scaling(Nterms))
297 2490 : energy_scaling = rvec !? just wondering, should this add up to 1, in which case we should cpp?
298 :
299 498 : CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", n_rep_val=n_rep)
300 498 : CPASSERT(n_rep == Nterms)
301 498 : CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=1, i_vals=ivec)
302 498 : Nelectron = SIZE(ivec)
303 498 : CPASSERT(Nelectron == k_alpha - k_beta)
304 1992 : ALLOCATE (occupations(2, Nelectron, Nterms))
305 7470 : occupations = 0
306 1494 : DO iterm = 1, Nterms
307 996 : CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=iterm, i_vals=ivec)
308 996 : CPASSERT(Nelectron == SIZE(ivec))
309 5976 : in_range = ALL(ivec >= 1) .AND. ALL(ivec <= 2)
310 996 : CPASSERT(in_range)
311 3486 : DO k = 1, Nelectron
312 2988 : occupations(ivec(k), k, iterm) = 1
313 : END DO
314 : END DO
315 :
316 : ! set up general data structures
317 : ! density matrices, kohn-sham matrices
318 :
319 498 : NULLIFY (matrix_p)
320 498 : CALL dbcsr_allocate_matrix_set(matrix_p, Nspin)
321 1494 : DO ispin = 1, Nspin
322 996 : ALLOCATE (matrix_p(ispin)%matrix)
323 : CALL dbcsr_copy(matrix_p(ispin)%matrix, rho_ao(1)%matrix, &
324 996 : name="density matrix low spin roks")
325 1494 : CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
326 : END DO
327 :
328 498 : NULLIFY (matrix_h)
329 498 : IF (.NOT. use_gauxc) THEN
330 466 : CALL dbcsr_allocate_matrix_set(matrix_h, Nspin)
331 1398 : DO ispin = 1, Nspin
332 932 : ALLOCATE (matrix_h(ispin)%matrix)
333 : CALL dbcsr_copy(matrix_h(ispin)%matrix, rho_ao(1)%matrix, &
334 932 : name="KS matrix low spin roks")
335 1398 : CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
336 : END DO
337 : END IF
338 :
339 498 : IF (do_hfx) THEN
340 250 : NULLIFY (matrix_hfx)
341 250 : CALL dbcsr_allocate_matrix_set(matrix_hfx, Nspin)
342 750 : DO ispin = 1, Nspin
343 500 : ALLOCATE (matrix_hfx(ispin)%matrix)
344 : CALL dbcsr_copy(matrix_hfx(ispin)%matrix, rho_ao(1)%matrix, &
345 750 : name="HFX matrix low spin roks")
346 : END DO
347 : END IF
348 :
349 : ! grids in real and g space for rho and vxc
350 : ! GAPW meta-GGAs need the kinetic-energy density of each configuration.
351 498 : NULLIFY (rho_r, rho_g, tau, tau_g, vxc_tau, vxc)
352 498 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
353 :
354 498 : IF (.NOT. use_gauxc) THEN
355 1398 : ALLOCATE (rho_r(Nspin))
356 1398 : ALLOCATE (rho_g(Nspin))
357 626 : IF (compute_tau) ALLOCATE (tau(Nspin), tau_g(Nspin))
358 1398 : DO ispin = 1, Nspin
359 932 : CALL auxbas_pw_pool%create_pw(rho_r(ispin))
360 932 : CALL auxbas_pw_pool%create_pw(rho_g(ispin))
361 1398 : IF (compute_tau) THEN
362 80 : CALL auxbas_pw_pool%create_pw(tau(ispin))
363 80 : CALL auxbas_pw_pool%create_pw(tau_g(ispin))
364 : END IF
365 : END DO
366 466 : CALL auxbas_pw_pool%create_pw(work_v_rspace)
367 : END IF
368 :
369 : ! get mo matrices needed to construct the density matrices
370 : ! we will base all on the alpha spin matrix, obviously possible in ROKS
371 498 : CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
372 498 : NULLIFY (fm_scaled, fm_deriv)
373 498 : CALL dbcsr_init_p(fm_scaled)
374 498 : CALL dbcsr_init_p(fm_deriv)
375 498 : CALL dbcsr_copy(fm_scaled, mo_coeff)
376 498 : CALL dbcsr_copy(fm_deriv, mo_coeff)
377 :
378 1494 : ALLOCATE (scaling(k_alpha))
379 :
380 : ! for each term, add it with the given scaling factor to the energy, and compute the required derivatives
381 1494 : DO iterm = 1, Nterms
382 :
383 2988 : DO ispin = 1, Nspin
384 : ! compute the proper density matrices with the required occupations
385 1992 : CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
386 17024 : scaling = 1.0_dp
387 5976 : scaling(k_alpha - Nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
388 1992 : CALL dbcsr_copy(fm_scaled, mo_coeff)
389 1992 : CALL dbcsr_scale_by_vector(fm_scaled, scaling, side='right')
390 : CALL dbcsr_multiply('n', 't', 1.0_dp, mo_coeff, fm_scaled, &
391 1992 : 0.0_dp, matrix_p(ispin)%matrix, retain_sparsity=.TRUE.)
392 2988 : IF (.NOT. use_gauxc) THEN
393 : ! compute the densities on the grid
394 : CALL calculate_rho_elec(matrix_p=matrix_p(ispin)%matrix, &
395 : rho=rho_r(ispin), rho_gspace=rho_g(ispin), &
396 1864 : ks_env=ks_env, soft_valid=gapw)
397 1864 : IF (compute_tau) THEN
398 : CALL calculate_rho_elec(matrix_p=matrix_p(ispin)%matrix, &
399 : rho=tau(ispin), rho_gspace=tau_g(ispin), &
400 160 : ks_env=ks_env, soft_valid=gapw, compute_tau=.TRUE.)
401 : END IF
402 : END IF
403 : END DO
404 :
405 996 : IF (gapw .AND. .NOT. use_gauxc) THEN
406 208 : CALL low_spin_roks_gapw_xc(qs_env, matrix_p, rho_atom_set, xc_section, just_energy, exc1)
407 208 : energy%exc1 = energy%exc1 + energy_scaling(iterm)*exc1
408 : END IF
409 :
410 : ! GauXC evaluates each configuration on the complete molecular AO density.
411 : ! Its matrix already contains all rho/gradient/tau/model-chain derivatives.
412 932 : IF (use_gauxc) THEN
413 64 : IF (ASSOCIATED(matrix_h)) CALL dbcsr_deallocate_matrix_set(matrix_h)
414 64 : matrix_p2(1:Nspin, 1:1) => matrix_p(1:Nspin)
415 : CALL apply_gauxc(qs_env, xc_section, .FALSE., rho_ao_external=matrix_p2, &
416 64 : exc_external=exc, matrix_vxc_external=matrix_h)
417 932 : ELSE IF (just_energy) THEN
418 : exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
419 168 : weights=weights, pw_pool=xc_pw_pool)
420 : ELSE
421 764 : CPASSERT(.NOT. compute_virial)
422 : CALL xc_vxc_pw_create(vxc_rho=vxc, rho_r=rho_r, &
423 : rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
424 : weights=weights, pw_pool=xc_pw_pool, &
425 764 : compute_virial=.FALSE., virial_xc=virial_xc_tmp)
426 : END IF
427 :
428 996 : energy%exc = energy%exc + energy_scaling(iterm)*exc
429 :
430 996 : IF (do_hfx) THEN
431 : ! Add Hartree-Fock contribution
432 1500 : DO ispin = 1, Nspin
433 1500 : CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
434 : END DO
435 500 : ehfx = energy%ex
436 : CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, &
437 500 : recalc_integrals=.FALSE., update_energy=.TRUE.)
438 500 : energy%ex = ehfx + energy_scaling(iterm)*energy%ex
439 : END IF
440 :
441 : ! add the corresponding derivatives to the MO derivatives
442 1494 : IF (.NOT. just_energy) THEN
443 804 : IF (use_gauxc) THEN
444 120 : DO ispin = 1, Nspin
445 120 : CALL dbcsr_scale(matrix_h(ispin)%matrix, energy_scaling(iterm))
446 : END DO
447 : ELSE
448 : ! get the potential in matrix form
449 2292 : DO ispin = 1, Nspin
450 1528 : CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
451 : ! use a work_v_rspace
452 1528 : CALL pw_axpy(vxc(ispin), work_v_rspace, energy_scaling(iterm)*vxc(ispin)%pw_grid%dvol, 0.0_dp)
453 : CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), &
454 1528 : qs_env=qs_env, calculate_forces=calculate_forces, gapw=gapw)
455 1528 : CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
456 2292 : IF (ASSOCIATED(vxc_tau)) THEN
457 : CALL pw_axpy(vxc_tau(ispin), work_v_rspace, &
458 96 : energy_scaling(iterm)*vxc_tau(ispin)%pw_grid%dvol, 0.0_dp)
459 : CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), &
460 96 : qs_env=qs_env, calculate_forces=.FALSE., gapw=gapw, compute_tau=.TRUE.)
461 96 : CALL auxbas_pw_pool%give_back_pw(vxc_tau(ispin))
462 : END IF
463 : END DO
464 764 : DEALLOCATE (vxc)
465 764 : IF (ASSOCIATED(vxc_tau)) DEALLOCATE (vxc_tau)
466 :
467 764 : IF (gapw) THEN
468 : CALL update_ks_atom(qs_env, matrix_h, matrix_p, .FALSE., &
469 128 : rho_atom_external=rho_atom_set, kscale=energy_scaling(iterm))
470 : END IF
471 : END IF
472 :
473 804 : IF (do_hfx) THEN
474 : ! add HFX contribution
475 1212 : DO ispin = 1, Nspin
476 : CALL dbcsr_add(matrix_h(ispin)%matrix, matrix_hfx(ispin)%matrix, &
477 1212 : 1.0_dp, energy_scaling(iterm))
478 : END DO
479 404 : IF (calculate_forces) THEN
480 8 : CALL get_qs_env(qs_env, x_data=x_data, para_env=para_env)
481 8 : IF (x_data(1, 1)%n_rep_hf /= 1) THEN
482 : CALL cp_abort(__LOCATION__, "Multiple HFX section forces not compatible "// &
483 0 : "with low spin ROKS method.")
484 : END IF
485 8 : IF (x_data(1, 1)%do_hfx_ri) THEN
486 0 : CALL cp_abort(__LOCATION__, "HFX_RI forces not compatible with low spin ROKS method.")
487 : ELSE
488 8 : irep = 1
489 8 : NULLIFY (mdummy)
490 8 : matrix_p2(1:Nspin, 1:1) => matrix_p(1:Nspin)
491 : CALL derivatives_four_center(qs_env, matrix_p2, mdummy, hfx_section, para_env, &
492 : irep, compute_virial, &
493 8 : adiabatic_rescale_factor=energy_scaling(iterm))
494 : END IF
495 : END IF
496 :
497 : END IF
498 :
499 : ! add this to the mo_derivs, again based on the alpha mo_coeff
500 2412 : DO ispin = 1, Nspin
501 : CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, &
502 1608 : 0.0_dp, dbcsr_deriv, last_column=k_alpha)
503 :
504 13536 : scaling = 1.0_dp
505 4824 : scaling(k_alpha - Nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
506 1608 : CALL dbcsr_scale_by_vector(dbcsr_deriv, scaling, side='right')
507 2412 : CALL dbcsr_add(mo_derivs(1)%matrix, dbcsr_deriv, 1.0_dp, 1.0_dp)
508 : END DO
509 :
510 : END IF
511 :
512 : END DO
513 :
514 : ! release allocated memory
515 498 : IF (.NOT. use_gauxc) THEN
516 1398 : DO ispin = 1, Nspin
517 932 : CALL auxbas_pw_pool%give_back_pw(rho_r(ispin))
518 932 : CALL auxbas_pw_pool%give_back_pw(rho_g(ispin))
519 1398 : IF (compute_tau) THEN
520 80 : CALL auxbas_pw_pool%give_back_pw(tau(ispin))
521 80 : CALL auxbas_pw_pool%give_back_pw(tau_g(ispin))
522 : END IF
523 : END DO
524 466 : DEALLOCATE (rho_r, rho_g)
525 466 : IF (compute_tau) DEALLOCATE (tau, tau_g)
526 466 : CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
527 : END IF
528 498 : CALL dbcsr_deallocate_matrix_set(matrix_p)
529 498 : IF (ASSOCIATED(matrix_h)) CALL dbcsr_deallocate_matrix_set(matrix_h)
530 498 : IF (do_hfx) THEN
531 250 : CALL dbcsr_deallocate_matrix_set(matrix_hfx)
532 : END IF
533 :
534 498 : CALL dbcsr_release_p(fm_deriv)
535 498 : CALL dbcsr_release_p(fm_scaled)
536 :
537 498 : DEALLOCATE (occupations)
538 498 : DEALLOCATE (energy_scaling)
539 498 : DEALLOCATE (scaling)
540 :
541 498 : CALL dbcsr_release_p(dbcsr_deriv)
542 :
543 498 : IF (ASSOCIATED(rho_atom_set)) CALL deallocate_rho_atom_set(rho_atom_set)
544 :
545 498 : CALL timestop(handle)
546 :
547 138651 : END SUBROUTINE low_spin_roks
548 :
549 : ! **************************************************************************************************
550 : !> \brief Evaluate the atomic XC correction for one low-spin ROKS configuration.
551 : !> Keep atomic densities and potential integrals separate from the reference
552 : !> state. All configurations share the total density, so no electrostatic
553 : !> correction or compensation-charge update is needed here.
554 : !> \param qs_env Quickstep environment
555 : !> \param matrix_p Spin density matrices of the configuration
556 : !> \param rho_atom_set Scratch atomic densities and XC potential integrals
557 : !> \param xc_section Exchange-correlation input
558 : !> \param just_energy Skip potential integrals when only the energy is requested
559 : !> \param exc1 Atomic hard-minus-soft XC energy
560 : ! **************************************************************************************************
561 208 : SUBROUTINE low_spin_roks_gapw_xc(qs_env, matrix_p, rho_atom_set, xc_section, just_energy, exc1)
562 : TYPE(qs_environment_type), POINTER :: qs_env
563 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_p
564 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
565 : TYPE(section_vals_type), POINTER :: xc_section
566 : LOGICAL, INTENT(IN) :: just_energy
567 : REAL(KIND=dp), INTENT(OUT) :: exc1
568 :
569 : INTEGER :: ikind, natom
570 208 : INTEGER, DIMENSION(:), POINTER :: atom_list
571 : LOGICAL :: paw_atom
572 : REAL(KIND=dp) :: hard_abs_spin, hard_spin, soft_abs_spin, &
573 : soft_spin
574 : REAL(KIND=dp), DIMENSION(2) :: hard_total, soft_total
575 208 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
576 208 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p2
577 : TYPE(dft_control_type), POINTER :: dft_control
578 : TYPE(mp_para_env_type), POINTER :: para_env
579 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
580 208 : POINTER :: sab
581 : TYPE(oce_matrix_type), POINTER :: oce
582 208 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
583 :
584 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
585 208 : dft_control=dft_control, para_env=para_env, oce=oce, sab_orb=sab)
586 208 : CPASSERT(dft_control%nspins == 2)
587 208 : CPASSERT(dft_control%nimages == 1)
588 208 : IF (.NOT. ASSOCIATED(rho_atom_set)) THEN
589 104 : CALL allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
590 : END IF
591 208 : matrix_p2(1:2, 1:1) => matrix_p
592 208 : CALL calculate_rho_atom_coeff(qs_env, matrix_p2, rho_atom_set, qs_kind_set, oce, sab, para_env)
593 :
594 : ! Only local diagnostics: do not overwrite the reference charge summaries.
595 208 : hard_total = 0.0_dp
596 208 : soft_total = 0.0_dp
597 208 : hard_spin = 0.0_dp
598 208 : soft_spin = 0.0_dp
599 208 : hard_abs_spin = 0.0_dp
600 208 : soft_abs_spin = 0.0_dp
601 416 : DO ikind = 1, SIZE(atomic_kind_set)
602 208 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
603 208 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
604 416 : IF (paw_atom) THEN
605 : CALL calculate_rho_atom(para_env, rho_atom_set, qs_kind_set(ikind), atom_list, natom, 2, &
606 208 : hard_total, soft_total, hard_spin, soft_spin, hard_abs_spin, soft_abs_spin)
607 : END IF
608 : END DO
609 208 : CALL zero_rho_atom_integrals(rho_atom_set)
610 208 : exc1 = 0.0_dp
611 : CALL calculate_vxc_atom(qs_env, just_energy, exc1, rho_atom_set_external=rho_atom_set, &
612 208 : xc_section_external=xc_section, calculate_forces=.FALSE.)
613 208 : END SUBROUTINE low_spin_roks_gapw_xc
614 :
615 : ! **************************************************************************************************
616 : !> \brief do sic calculations on explicit orbitals
617 : !> \param energy ...
618 : !> \param qs_env ...
619 : !> \param dft_control ...
620 : !> \param poisson_env ...
621 : !> \param just_energy ...
622 : !> \param calculate_forces ...
623 : !> \param auxbas_pw_pool ...
624 : ! **************************************************************************************************
625 136161 : SUBROUTINE sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, &
626 : calculate_forces, auxbas_pw_pool)
627 :
628 : TYPE(qs_energy_type), POINTER :: energy
629 : TYPE(qs_environment_type), POINTER :: qs_env
630 : TYPE(dft_control_type), POINTER :: dft_control
631 : TYPE(pw_poisson_type), POINTER :: poisson_env
632 : LOGICAL, INTENT(IN) :: just_energy, calculate_forces
633 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
634 :
635 : CHARACTER(*), PARAMETER :: routineN = 'sic_explicit_orbitals'
636 :
637 : INTEGER :: handle, i, Iorb, k_alpha, k_beta, Norb
638 136161 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: sic_orbital_list
639 : LOGICAL :: compute_virial, uniform_occupation
640 : REAL(KIND=dp) :: ener, exc
641 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp
642 : TYPE(cell_type), POINTER :: cell
643 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
644 : TYPE(cp_fm_type) :: matrix_hv, matrix_v
645 136161 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_derivs_local
646 : TYPE(cp_fm_type), POINTER :: mo_coeff
647 : TYPE(dbcsr_p_type) :: orb_density_matrix_p, orb_h_p
648 136161 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, rho_ao, tmp_dbcsr
649 : TYPE(dbcsr_type), POINTER :: orb_density_matrix, orb_h
650 136161 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
651 : TYPE(pw_c1d_gs_type) :: work_v_gspace
652 136161 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
653 : TYPE(pw_c1d_gs_type), TARGET :: orb_rho_g, tmp_g
654 : TYPE(pw_env_type), POINTER :: pw_env
655 : TYPE(pw_pool_type), POINTER :: xc_pw_pool
656 : TYPE(pw_r3d_rs_type) :: work_v_rspace
657 136161 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc, vxc_tau
658 : TYPE(pw_r3d_rs_type), POINTER :: weights
659 : TYPE(pw_r3d_rs_type), TARGET :: orb_rho_r, tmp_r
660 : TYPE(qs_ks_env_type), POINTER :: ks_env
661 : TYPE(qs_rho_type), POINTER :: rho
662 : TYPE(section_vals_type), POINTER :: input, xc_section
663 : TYPE(virial_type), POINTER :: virial
664 :
665 136161 : IF (dft_control%sic_method_id /= sic_eo) RETURN
666 :
667 0 : CALL timeset(routineN, handle)
668 :
669 0 : NULLIFY (tau, vxc_tau, mo_derivs, ks_env, rho_ao)
670 :
671 : ! generate the lists of orbitals that need sic treatment
672 : CALL get_qs_env(qs_env, &
673 : ks_env=ks_env, &
674 : mo_derivs=mo_derivs, &
675 : mos=mo_array, &
676 : rho=rho, &
677 : xcint_weights=weights, &
678 : pw_env=pw_env, &
679 : input=input, &
680 : cell=cell, &
681 0 : virial=virial)
682 :
683 0 : CALL qs_rho_get(rho, rho_ao=rho_ao)
684 :
685 0 : compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
686 0 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
687 :
688 0 : DO i = 1, SIZE(mo_array) !fm->dbcsr
689 0 : IF (mo_array(i)%use_mo_coeff_b) THEN !fm->dbcsr
690 : CALL copy_dbcsr_to_fm(mo_array(i)%mo_coeff_b, &
691 0 : mo_array(i)%mo_coeff) !fm->dbcsr
692 : END IF !fm->dbcsr
693 : END DO !fm->dbcsr
694 :
695 0 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
696 :
697 : ! we have two spins
698 0 : CPASSERT(SIZE(mo_array, 1) == 2)
699 : ! we want uniform occupations
700 0 : CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
701 0 : CPASSERT(uniform_occupation)
702 0 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff, uniform_occupation=uniform_occupation)
703 0 : CPASSERT(uniform_occupation)
704 :
705 0 : NULLIFY (tmp_dbcsr)
706 0 : CALL dbcsr_allocate_matrix_set(tmp_dbcsr, SIZE(mo_derivs, 1))
707 0 : DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
708 : !
709 0 : NULLIFY (tmp_dbcsr(i)%matrix)
710 0 : CALL dbcsr_init_p(tmp_dbcsr(i)%matrix)
711 0 : CALL dbcsr_copy(tmp_dbcsr(i)%matrix, mo_derivs(i)%matrix)
712 0 : CALL dbcsr_set(tmp_dbcsr(i)%matrix, 0.0_dp)
713 : END DO !fm->dbcsr
714 :
715 0 : k_alpha = 0; k_beta = 0
716 0 : SELECT CASE (dft_control%sic_list_id)
717 : CASE (sic_list_all)
718 :
719 0 : CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
720 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
721 :
722 0 : IF (SIZE(mo_array, 1) > 1) THEN
723 0 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
724 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
725 : END IF
726 :
727 0 : Norb = k_alpha + k_beta
728 0 : ALLOCATE (sic_orbital_list(3, Norb))
729 :
730 0 : iorb = 0
731 0 : DO i = 1, k_alpha
732 0 : iorb = iorb + 1
733 0 : sic_orbital_list(1, iorb) = 1
734 0 : sic_orbital_list(2, iorb) = i
735 0 : sic_orbital_list(3, iorb) = 1
736 : END DO
737 0 : DO i = 1, k_beta
738 0 : iorb = iorb + 1
739 0 : sic_orbital_list(1, iorb) = 2
740 0 : sic_orbital_list(2, iorb) = i
741 0 : IF (SIZE(mo_derivs, 1) == 1) THEN
742 0 : sic_orbital_list(3, iorb) = 1
743 : ELSE
744 0 : sic_orbital_list(3, iorb) = 2
745 : END IF
746 : END DO
747 :
748 : CASE (sic_list_unpaired)
749 : ! we have two spins
750 0 : CPASSERT(SIZE(mo_array, 1) == 2)
751 : ! we have them restricted
752 0 : CPASSERT(SIZE(mo_derivs, 1) == 1)
753 0 : CPASSERT(dft_control%restricted)
754 :
755 0 : CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
756 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
757 :
758 0 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
759 0 : CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
760 :
761 0 : Norb = k_alpha - k_beta
762 0 : ALLOCATE (sic_orbital_list(3, Norb))
763 :
764 0 : iorb = 0
765 0 : DO i = k_beta + 1, k_alpha
766 0 : iorb = iorb + 1
767 0 : sic_orbital_list(1, iorb) = 1
768 0 : sic_orbital_list(2, iorb) = i
769 : ! we are guaranteed to be restricted
770 0 : sic_orbital_list(3, iorb) = 1
771 : END DO
772 :
773 : CASE DEFAULT
774 0 : CPABORT("Unknown dft_control%sic_list_id")
775 : END SELECT
776 :
777 : ! data needed for each of the orbs
778 0 : CALL auxbas_pw_pool%create_pw(orb_rho_r)
779 0 : CALL auxbas_pw_pool%create_pw(tmp_r)
780 0 : CALL auxbas_pw_pool%create_pw(orb_rho_g)
781 0 : CALL auxbas_pw_pool%create_pw(tmp_g)
782 0 : CALL auxbas_pw_pool%create_pw(work_v_gspace)
783 0 : CALL auxbas_pw_pool%create_pw(work_v_rspace)
784 :
785 0 : ALLOCATE (orb_density_matrix)
786 : CALL dbcsr_copy(orb_density_matrix, rho_ao(1)%matrix, &
787 0 : name="orb_density_matrix")
788 0 : CALL dbcsr_set(orb_density_matrix, 0.0_dp)
789 0 : orb_density_matrix_p%matrix => orb_density_matrix
790 :
791 0 : ALLOCATE (orb_h)
792 : CALL dbcsr_copy(orb_h, rho_ao(1)%matrix, &
793 0 : name="orb_density_matrix")
794 0 : CALL dbcsr_set(orb_h, 0.0_dp)
795 0 : orb_h_p%matrix => orb_h
796 :
797 0 : CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
798 :
799 : CALL cp_fm_struct_create(fm_struct_tmp, ncol_global=1, &
800 0 : template_fmstruct=mo_coeff%matrix_struct)
801 0 : CALL cp_fm_create(matrix_v, fm_struct_tmp, name="matrix_v")
802 0 : CALL cp_fm_create(matrix_hv, fm_struct_tmp, name="matrix_hv")
803 0 : CALL cp_fm_struct_release(fm_struct_tmp)
804 :
805 0 : ALLOCATE (mo_derivs_local(SIZE(mo_array, 1)))
806 0 : DO I = 1, SIZE(mo_array, 1)
807 0 : CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff)
808 0 : CALL cp_fm_create(mo_derivs_local(I), mo_coeff%matrix_struct)
809 : END DO
810 :
811 0 : ALLOCATE (rho_r(2))
812 0 : rho_r(1) = orb_rho_r
813 0 : rho_r(2) = tmp_r
814 0 : CALL pw_zero(tmp_r)
815 :
816 0 : ALLOCATE (rho_g(2))
817 0 : rho_g(1) = orb_rho_g
818 0 : rho_g(2) = tmp_g
819 0 : CALL pw_zero(tmp_g)
820 :
821 0 : NULLIFY (vxc)
822 : ! now apply to SIC correction to each selected orbital
823 0 : DO iorb = 1, Norb
824 : ! extract the proper orbital from the mo_coeff
825 0 : CALL get_mo_set(mo_set=mo_array(sic_orbital_list(1, iorb)), mo_coeff=mo_coeff)
826 0 : CALL cp_fm_to_fm(mo_coeff, matrix_v, 1, sic_orbital_list(2, iorb), 1)
827 :
828 : ! construct the density matrix and the corresponding density
829 0 : CALL dbcsr_set(orb_density_matrix, 0.0_dp)
830 : CALL cp_dbcsr_plus_fm_fm_t(orb_density_matrix, matrix_v=matrix_v, ncol=1, &
831 0 : alpha=1.0_dp)
832 :
833 : CALL calculate_rho_elec(matrix_p=orb_density_matrix, &
834 : rho=orb_rho_r, rho_gspace=orb_rho_g, &
835 0 : ks_env=ks_env)
836 :
837 : ! compute the energy functional for this orbital and its derivative
838 :
839 0 : CALL pw_poisson_solve(poisson_env, orb_rho_g, ener, work_v_gspace)
840 : ! no PBC correction is done here, see "calc_v_sic_rspace" for SIC methods
841 : ! with PBC aware corrections
842 0 : energy%hartree = energy%hartree - dft_control%sic_scaling_a*ener
843 0 : IF (.NOT. just_energy) THEN
844 0 : CALL pw_transfer(work_v_gspace, work_v_rspace)
845 0 : CALL pw_scale(work_v_rspace, -dft_control%sic_scaling_a*work_v_rspace%pw_grid%dvol)
846 0 : CALL dbcsr_set(orb_h, 0.0_dp)
847 : END IF
848 :
849 0 : IF (just_energy) THEN
850 : exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
851 0 : weights=weights, pw_pool=xc_pw_pool)
852 : ELSE
853 0 : CPASSERT(.NOT. compute_virial)
854 : CALL xc_vxc_pw_create(vxc_rho=vxc, rho_r=rho_r, &
855 : rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
856 : weights=weights, pw_pool=xc_pw_pool, &
857 0 : compute_virial=compute_virial, virial_xc=virial_xc_tmp)
858 : ! add to the existing work_v_rspace
859 0 : CALL pw_axpy(vxc(1), work_v_rspace, -dft_control%sic_scaling_b*vxc(1)%pw_grid%dvol)
860 : END IF
861 0 : energy%exc = energy%exc - dft_control%sic_scaling_b*exc
862 :
863 0 : IF (.NOT. just_energy) THEN
864 : ! note, orb_h (which is being pointed to with orb_h_p) is zeroed above
865 : CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=orb_density_matrix_p, hmat=orb_h_p, &
866 0 : qs_env=qs_env, calculate_forces=calculate_forces)
867 :
868 : ! add this to the mo_derivs
869 0 : CALL cp_dbcsr_sm_fm_multiply(orb_h, matrix_v, matrix_hv, 1)
870 : ! silly trick, copy to an array of the right size and add to mo_derivs
871 0 : CALL cp_fm_set_all(mo_derivs_local(sic_orbital_list(3, iorb)), 0.0_dp)
872 0 : CALL cp_fm_to_fm(matrix_hv, mo_derivs_local(sic_orbital_list(3, iorb)), 1, 1, sic_orbital_list(2, iorb))
873 : CALL copy_fm_to_dbcsr(mo_derivs_local(sic_orbital_list(3, iorb)), &
874 0 : tmp_dbcsr(sic_orbital_list(3, iorb))%matrix)
875 : CALL dbcsr_add(mo_derivs(sic_orbital_list(3, iorb))%matrix, &
876 0 : tmp_dbcsr(sic_orbital_list(3, iorb))%matrix, 1.0_dp, 1.0_dp)
877 : !
878 : ! need to deallocate vxc
879 0 : CALL xc_pw_pool%give_back_pw(vxc(1))
880 0 : CALL xc_pw_pool%give_back_pw(vxc(2))
881 0 : DEALLOCATE (vxc)
882 :
883 : END IF
884 :
885 : END DO
886 :
887 0 : CALL auxbas_pw_pool%give_back_pw(orb_rho_r)
888 0 : CALL auxbas_pw_pool%give_back_pw(tmp_r)
889 0 : CALL auxbas_pw_pool%give_back_pw(orb_rho_g)
890 0 : CALL auxbas_pw_pool%give_back_pw(tmp_g)
891 0 : CALL auxbas_pw_pool%give_back_pw(work_v_gspace)
892 0 : CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
893 :
894 0 : CALL dbcsr_deallocate_matrix(orb_density_matrix)
895 0 : CALL dbcsr_deallocate_matrix(orb_h)
896 0 : CALL cp_fm_release(matrix_v)
897 0 : CALL cp_fm_release(matrix_hv)
898 0 : CALL cp_fm_release(mo_derivs_local)
899 0 : DEALLOCATE (rho_r)
900 0 : DEALLOCATE (rho_g)
901 :
902 0 : CALL dbcsr_deallocate_matrix_set(tmp_dbcsr) !fm->dbcsr
903 :
904 0 : CALL timestop(handle)
905 :
906 136161 : END SUBROUTINE sic_explicit_orbitals
907 :
908 : ! **************************************************************************************************
909 : !> \brief do sic calculations on the spin density
910 : !> \param v_sic_rspace ...
911 : !> \param energy ...
912 : !> \param qs_env ...
913 : !> \param dft_control ...
914 : !> \param rho ...
915 : !> \param poisson_env ...
916 : !> \param just_energy ...
917 : !> \param calculate_forces ...
918 : !> \param auxbas_pw_pool ...
919 : ! **************************************************************************************************
920 136161 : SUBROUTINE calc_v_sic_rspace(v_sic_rspace, energy, &
921 : qs_env, dft_control, rho, poisson_env, just_energy, &
922 : calculate_forces, auxbas_pw_pool)
923 :
924 : TYPE(pw_r3d_rs_type), POINTER :: v_sic_rspace
925 : TYPE(qs_energy_type), POINTER :: energy
926 : TYPE(qs_environment_type), POINTER :: qs_env
927 : TYPE(dft_control_type), POINTER :: dft_control
928 : TYPE(qs_rho_type), POINTER :: rho
929 : TYPE(pw_poisson_type), POINTER :: poisson_env
930 : LOGICAL, INTENT(IN) :: just_energy, calculate_forces
931 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
932 :
933 : INTEGER :: i, nelec, nelec_a, nelec_b, nforce
934 : REAL(kind=dp) :: ener, full_scaling, scaling
935 136161 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: store_forces
936 136161 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
937 : TYPE(pw_c1d_gs_type) :: work_rho, work_v
938 136161 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
939 136161 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
940 :
941 136161 : NULLIFY (mo_array, rho_g)
942 :
943 136161 : IF (dft_control%sic_method_id == sic_none) RETURN
944 94 : IF (dft_control%sic_method_id == sic_eo) RETURN
945 :
946 94 : IF (dft_control%qs_control%gapw) THEN
947 0 : CPABORT("sic and GAPW not yet compatible")
948 : END IF
949 :
950 : ! OK, right now we like two spins to do sic, could be relaxed for AD
951 94 : CPASSERT(dft_control%nspins == 2)
952 :
953 94 : CALL auxbas_pw_pool%create_pw(work_rho)
954 94 : CALL auxbas_pw_pool%create_pw(work_v)
955 :
956 94 : CALL qs_rho_get(rho, rho_g=rho_g)
957 :
958 : ! Hartree sic corrections
959 188 : SELECT CASE (dft_control%sic_method_id)
960 : CASE (sic_mauri_us, sic_mauri_spz)
961 94 : CALL pw_copy(rho_g(1), work_rho)
962 94 : CALL pw_axpy(rho_g(2), work_rho, alpha=-1._dp)
963 94 : CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
964 : CASE (sic_ad)
965 : ! find out how many elecs we have
966 0 : CALL get_qs_env(qs_env, mos=mo_array)
967 0 : CALL get_mo_set(mo_set=mo_array(1), nelectron=nelec_a)
968 0 : CALL get_mo_set(mo_set=mo_array(2), nelectron=nelec_b)
969 0 : nelec = nelec_a + nelec_b
970 0 : CALL pw_copy(rho_g(1), work_rho)
971 0 : CALL pw_axpy(rho_g(2), work_rho)
972 0 : scaling = 1.0_dp/REAL(nelec, KIND=dp)
973 0 : CALL pw_scale(work_rho, scaling)
974 0 : CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
975 : CASE DEFAULT
976 188 : CPABORT("Unknown sic method id")
977 : END SELECT
978 :
979 : ! Correct for DDAP charges (if any)
980 : ! storing whatever force might be there from previous decoupling
981 94 : IF (calculate_forces) THEN
982 12 : CALL get_qs_env(qs_env=qs_env, force=force)
983 12 : nforce = 0
984 36 : DO i = 1, SIZE(force)
985 36 : nforce = nforce + SIZE(force(i)%ch_pulay, 2)
986 : END DO
987 36 : ALLOCATE (store_forces(3, nforce))
988 36 : nforce = 0
989 36 : DO i = 1, SIZE(force)
990 312 : store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2)) = force(i)%ch_pulay(:, :)
991 312 : force(i)%ch_pulay(:, :) = 0.0_dp
992 36 : nforce = nforce + SIZE(force(i)%ch_pulay, 2)
993 : END DO
994 : END IF
995 :
996 : CALL cp_ddapc_apply_CD(qs_env, &
997 : work_rho, &
998 : ener, &
999 : v_hartree_gspace=work_v, &
1000 : calculate_forces=calculate_forces, &
1001 94 : Itype_of_density="SPIN")
1002 :
1003 188 : SELECT CASE (dft_control%sic_method_id)
1004 : CASE (sic_mauri_us, sic_mauri_spz)
1005 94 : full_scaling = -dft_control%sic_scaling_a
1006 : CASE (sic_ad)
1007 0 : full_scaling = -dft_control%sic_scaling_a*nelec
1008 : CASE DEFAULT
1009 94 : CPABORT("Unknown sic method id")
1010 : END SELECT
1011 94 : energy%hartree = energy%hartree + full_scaling*ener
1012 :
1013 : ! add scaled forces, restoring the old
1014 94 : IF (calculate_forces) THEN
1015 12 : nforce = 0
1016 36 : DO i = 1, SIZE(force)
1017 : force(i)%ch_pulay(:, :) = force(i)%ch_pulay(:, :)*full_scaling + &
1018 312 : store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2))
1019 36 : nforce = nforce + SIZE(force(i)%ch_pulay, 2)
1020 : END DO
1021 : END IF
1022 :
1023 94 : IF (.NOT. just_energy) THEN
1024 60 : ALLOCATE (v_sic_rspace)
1025 60 : CALL auxbas_pw_pool%create_pw(v_sic_rspace)
1026 60 : CALL pw_transfer(work_v, v_sic_rspace)
1027 : ! also take into account the scaling (in addition to the volume element)
1028 : CALL pw_scale(v_sic_rspace, &
1029 60 : dft_control%sic_scaling_a*v_sic_rspace%pw_grid%dvol)
1030 : END IF
1031 :
1032 94 : CALL auxbas_pw_pool%give_back_pw(work_rho)
1033 94 : CALL auxbas_pw_pool%give_back_pw(work_v)
1034 :
1035 136255 : END SUBROUTINE calc_v_sic_rspace
1036 :
1037 : ! **************************************************************************************************
1038 : !> \brief Check each term of the energy and sum to total
1039 : !> \param energy ...
1040 : ! **************************************************************************************************
1041 136095 : SUBROUTINE check_sum_energies(energy)
1042 : TYPE(qs_energy_type), INTENT(INOUT), POINTER :: energy
1043 :
1044 : ! Reset the total energy counter, just in case
1045 : ! it is used to store some intermediate value
1046 136095 : energy%total = 0.0_dp
1047 :
1048 136095 : CALL check_sum_energy(energy%total, energy%core_overlap, "Core overlap energy")
1049 136095 : CALL check_sum_energy(energy%total, energy%core_self, "Core self energy")
1050 136095 : CALL check_sum_energy(energy%total, energy%core_cneo, "Quantum nuclear core energy")
1051 136095 : CALL check_sum_energy(energy%total, energy%core, "Core Hamiltonian energy")
1052 136095 : CALL check_sum_energy(energy%total, energy%hartree, "Hartree energy")
1053 136095 : CALL check_sum_energy(energy%total, energy%hartree_1c, "Energy from GAPW local Eh = 1 center integrals")
1054 136095 : CALL check_sum_energy(energy%total, energy%exc, "Exchange-correlation energy")
1055 136095 : CALL check_sum_energy(energy%total, energy%exc1, "Exc1 energy")
1056 136095 : CALL check_sum_energy(energy%total, energy%ex, "Ex energy")
1057 136095 : CALL check_sum_energy(energy%total, energy%dispersion, "Dispersion energy")
1058 136095 : CALL check_sum_energy(energy%total, energy%gcp, "gCP energy")
1059 136095 : CALL check_sum_energy(energy%total, energy%qmmm_el, "QM/MM Electrostatic energy")
1060 136095 : CALL check_sum_energy(energy%total, energy%mulliken, "Mulliken restraint energy")
1061 272294 : CALL check_sum_energy(energy%total, SUM(energy%ddapc_restraint), "DDAPC restraint energy")
1062 136095 : CALL check_sum_energy(energy%total, energy%s2_restraint, "S2 restraint energy")
1063 136095 : CALL check_sum_energy(energy%total, energy%dft_plus_u, "DFT+U energy")
1064 136095 : CALL check_sum_energy(energy%total, energy%kTS, "Electronic entropic contribution energy")
1065 136095 : CALL check_sum_energy(energy%total, energy%efield, "Electric field interaction energy")
1066 136095 : CALL check_sum_energy(energy%total, energy%efield_core, "Electric field core interaction energy")
1067 136095 : CALL check_sum_energy(energy%total, energy%ee, "External potential interaction energy")
1068 136095 : CALL check_sum_energy(energy%total, energy%ee_core, "External potential core interaction energy")
1069 136095 : CALL check_sum_energy(energy%total, energy%exc_aux_fit, "Wfn fit exchange-correlation energy")
1070 136095 : CALL check_sum_energy(energy%total, energy%image_charge, "QM/MM image charge energy energy")
1071 136095 : CALL check_sum_energy(energy%total, energy%sccs_pol, "SCCS polarisation energy")
1072 136095 : CALL check_sum_energy(energy%total, energy%cdft, "CDFT constraint energy")
1073 136095 : CALL check_sum_energy(energy%total, energy%exc1_aux_fit, "Wfn fit soft/hard atomic rho1 Exc contribution energy")
1074 136095 : CALL check_sum_energy(energy%total, energy%embed_corr, "Embedding potential correction energy")
1075 :
1076 136095 : IF (abnormal_value(energy%total)) THEN
1077 : CALL cp_abort(__LOCATION__, &
1078 0 : "Total energy after summing up the terms is an abnormal value (NaN/Inf).")
1079 : END IF
1080 :
1081 136095 : END SUBROUTINE check_sum_energies
1082 :
1083 : ! **************************************************************************************************
1084 : !> \brief Check one term of the energy and sum to total if normal
1085 : !> \param total ...
1086 : !> \param term ...
1087 : !> \param name ...
1088 : ! **************************************************************************************************
1089 3674565 : SUBROUTINE check_sum_energy(total, term, name)
1090 : REAL(KIND=dp), INTENT(INOUT) :: total
1091 : REAL(KIND=dp), INTENT(IN) :: term
1092 : CHARACTER(LEN=*), INTENT(IN) :: name
1093 :
1094 3674565 : IF (abnormal_value(term)) THEN
1095 : CALL cp_abort(__LOCATION__, &
1096 0 : TRIM(name)//" is an abnormal value (NaN/Inf).")
1097 : ELSE
1098 3674565 : total = total + term
1099 : END IF
1100 :
1101 3674565 : END SUBROUTINE check_sum_energy
1102 :
1103 : ! **************************************************************************************************
1104 : !> \brief ...
1105 : !> \param qs_env ...
1106 : !> \param rho ...
1107 : ! **************************************************************************************************
1108 272210 : SUBROUTINE print_densities(qs_env, rho)
1109 : TYPE(qs_environment_type), POINTER :: qs_env
1110 : TYPE(qs_rho_type), POINTER :: rho
1111 :
1112 : INTEGER :: img, ispin, n_electrons, output_unit
1113 : REAL(dp) :: tot1_h, tot1_s, tot_rho_r, trace, &
1114 : trace_tmp
1115 136105 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r_arr
1116 : TYPE(cell_type), POINTER :: cell
1117 : TYPE(cp_logger_type), POINTER :: logger
1118 136105 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, rho_ao
1119 : TYPE(dft_control_type), POINTER :: dft_control
1120 : TYPE(qs_charges_type), POINTER :: qs_charges
1121 136105 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1122 : TYPE(section_vals_type), POINTER :: input, scf_section
1123 :
1124 136105 : NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, &
1125 136105 : dft_control, tot_rho_r_arr, rho_ao)
1126 :
1127 272210 : logger => cp_get_default_logger()
1128 :
1129 : CALL get_qs_env(qs_env, &
1130 : qs_kind_set=qs_kind_set, &
1131 : cell=cell, qs_charges=qs_charges, &
1132 : input=input, &
1133 : matrix_s_kp=matrix_s, &
1134 136105 : dft_control=dft_control)
1135 :
1136 136105 : CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
1137 :
1138 136105 : scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
1139 : output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", &
1140 136105 : extension=".scfLog")
1141 :
1142 136105 : CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao)
1143 136105 : n_electrons = n_electrons - dft_control%charge
1144 136105 : tot_rho_r = accurate_sum(tot_rho_r_arr)
1145 :
1146 136105 : trace = 0
1147 136105 : IF (BTEST(cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES"), cp_p_file)) THEN
1148 4154 : DO ispin = 1, dft_control%nspins
1149 7526 : DO img = 1, dft_control%nimages
1150 3372 : CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp)
1151 5848 : trace = trace + trace_tmp
1152 : END DO
1153 : END DO
1154 : END IF
1155 :
1156 136105 : IF (output_unit > 0) THEN
1157 839 : WRITE (UNIT=output_unit, FMT="(/,T3,A,T41,F20.10)") "Trace(PS):", trace
1158 : WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
1159 839 : "Electronic density on regular grids: ", &
1160 839 : tot_rho_r, &
1161 : tot_rho_r + &
1162 839 : REAL(n_electrons, dp), &
1163 839 : "Core density on regular grids:", &
1164 839 : qs_charges%total_rho_core_rspace, &
1165 : qs_charges%total_rho_core_rspace + &
1166 : qs_charges%total_rho1_hard_nuc - &
1167 1678 : REAL(n_electrons + dft_control%charge, dp)
1168 : END IF
1169 136105 : IF (dft_control%qs_control%gapw) THEN
1170 23922 : tot1_h = qs_charges%total_rho1_hard(1)
1171 23922 : tot1_s = qs_charges%total_rho1_soft(1)
1172 28278 : DO ispin = 2, dft_control%nspins
1173 4356 : tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
1174 28278 : tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
1175 : END DO
1176 23922 : IF (output_unit > 0) THEN
1177 : WRITE (UNIT=output_unit, FMT="((T3,A,T41,2F20.10))") &
1178 399 : "Hard and soft densities (Lebedev):", &
1179 798 : tot1_h, tot1_s
1180 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
1181 399 : "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
1182 399 : tot_rho_r + tot1_h - tot1_s, &
1183 399 : "Total charge density (r-space): ", &
1184 : tot_rho_r + tot1_h - tot1_s &
1185 : + qs_charges%total_rho_core_rspace &
1186 798 : + qs_charges%total_rho1_hard_nuc
1187 399 : IF (qs_charges%total_rho1_hard_nuc /= 0.0_dp) THEN
1188 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
1189 0 : "Total CNEO nuc. char. den. (Lebedev): ", &
1190 0 : qs_charges%total_rho1_hard_nuc, &
1191 0 : "Total CNEO soft char. den. (Lebedev): ", &
1192 0 : qs_charges%total_rho1_soft_nuc_lebedev, &
1193 0 : "Total CNEO soft char. den. (r-space): ", &
1194 0 : qs_charges%total_rho1_soft_nuc_rspace, &
1195 0 : "Total soft Rho_e+n+0 (g-space):", &
1196 0 : qs_charges%total_rho_gspace
1197 : ELSE
1198 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
1199 399 : "Total Rho_soft + Rho0_soft (g-space):", &
1200 798 : qs_charges%total_rho_gspace
1201 : END IF
1202 : END IF
1203 : qs_charges%background = tot_rho_r + tot1_h - tot1_s + &
1204 : qs_charges%total_rho_core_rspace + &
1205 23922 : qs_charges%total_rho1_hard_nuc
1206 : ! only add total_rho1_hard_nuc for gapw as cneo requires gapw
1207 112183 : ELSE IF (dft_control%qs_control%gapw_xc) THEN
1208 4448 : tot1_h = qs_charges%total_rho1_hard(1)
1209 4448 : tot1_s = qs_charges%total_rho1_soft(1)
1210 4826 : DO ispin = 2, dft_control%nspins
1211 378 : tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
1212 4826 : tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
1213 : END DO
1214 4448 : IF (output_unit > 0) THEN
1215 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T41,2F20.10))") &
1216 0 : "Hard and soft densities (Lebedev):", &
1217 0 : tot1_h, tot1_s
1218 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
1219 0 : "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
1220 0 : accurate_sum(tot_rho_r_arr) + tot1_h - tot1_s
1221 : END IF
1222 : qs_charges%background = tot_rho_r + &
1223 4448 : qs_charges%total_rho_core_rspace
1224 : ELSE
1225 107735 : IF (output_unit > 0) THEN
1226 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
1227 440 : "Total charge density on r-space grids: ", &
1228 : tot_rho_r + &
1229 440 : qs_charges%total_rho_core_rspace, &
1230 440 : "Total charge density g-space grids: ", &
1231 880 : qs_charges%total_rho_gspace
1232 : END IF
1233 : qs_charges%background = tot_rho_r + &
1234 107735 : qs_charges%total_rho_core_rspace
1235 : END IF
1236 136105 : IF (output_unit > 0) WRITE (UNIT=output_unit, FMT="()")
1237 136105 : qs_charges%background = qs_charges%background/cell%deth
1238 :
1239 : CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1240 136105 : "PRINT%TOTAL_DENSITIES")
1241 :
1242 136105 : END SUBROUTINE print_densities
1243 :
1244 : ! **************************************************************************************************
1245 : !> \brief Print detailed energies
1246 : !>
1247 : !> \param qs_env ...
1248 : !> \param dft_control ...
1249 : !> \param input ...
1250 : !> \param energy ...
1251 : !> \param mulliken_order_p ...
1252 : !> \par History
1253 : !> refactoring 04.03.2011 [MI]
1254 : !> \author
1255 : ! **************************************************************************************************
1256 136105 : SUBROUTINE print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
1257 :
1258 : TYPE(qs_environment_type), POINTER :: qs_env
1259 : TYPE(dft_control_type), POINTER :: dft_control
1260 : TYPE(section_vals_type), POINTER :: input
1261 : TYPE(qs_energy_type), POINTER :: energy
1262 : REAL(KIND=dp), INTENT(IN) :: mulliken_order_p
1263 :
1264 : INTEGER :: bc, n, output_unit, psolver
1265 : REAL(KIND=dp) :: ddapc_order_p, implicit_ps_ehartree, &
1266 : s2_order_p
1267 : TYPE(cp_logger_type), POINTER :: logger
1268 : TYPE(pw_env_type), POINTER :: pw_env
1269 :
1270 136105 : logger => cp_get_default_logger()
1271 :
1272 136105 : NULLIFY (pw_env)
1273 136105 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1274 136105 : psolver = pw_env%poisson_env%parameters%solver
1275 :
1276 : output_unit = cp_print_key_unit_nr(logger, input, "DFT%SCF%PRINT%DETAILED_ENERGY", &
1277 136105 : extension=".scfLog")
1278 136105 : IF (output_unit > 0) THEN
1279 492 : IF (dft_control%do_admm) THEN
1280 : WRITE (UNIT=output_unit, FMT="((T3,A,T60,F20.10))") &
1281 0 : "Wfn fit exchange-correlation energy: ", energy%exc_aux_fit
1282 0 : IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1283 : WRITE (UNIT=output_unit, FMT="((T3,A,T60,F20.10))") &
1284 0 : "Wfn fit soft/hard atomic rho1 Exc contribution: ", energy%exc1_aux_fit
1285 : END IF
1286 : END IF
1287 492 : IF (dft_control%do_admm) THEN
1288 0 : IF (psolver == pw_poisson_implicit) THEN
1289 0 : implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1290 0 : bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1291 0 : SELECT CASE (bc)
1292 : CASE (MIXED_PERIODIC_BC, MIXED_BC)
1293 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1294 0 : "Core Hamiltonian energy: ", energy%core, &
1295 0 : "Hartree energy: ", implicit_ps_ehartree, &
1296 0 : "Electric enthalpy: ", energy%hartree, &
1297 0 : "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1298 : CASE (PERIODIC_BC, NEUMANN_BC)
1299 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1300 0 : "Core Hamiltonian energy: ", energy%core, &
1301 0 : "Hartree energy: ", energy%hartree, &
1302 0 : "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1303 : END SELECT
1304 : ELSE
1305 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1306 0 : "Core Hamiltonian energy: ", energy%core, &
1307 0 : "Hartree energy: ", energy%hartree, &
1308 0 : "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1309 : END IF
1310 : ELSE
1311 : !ZMP to print some variables at each step
1312 492 : IF (dft_control%apply_external_density) THEN
1313 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1314 0 : "DOING ZMP CALCULATION FROM EXTERNAL DENSITY "
1315 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1316 0 : "Core Hamiltonian energy: ", energy%core, &
1317 0 : "Hartree energy: ", energy%hartree
1318 492 : ELSE IF (dft_control%apply_external_vxc) THEN
1319 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1320 0 : "DOING ZMP READING EXTERNAL VXC "
1321 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1322 0 : "Core Hamiltonian energy: ", energy%core, &
1323 0 : "Hartree energy: ", energy%hartree
1324 : ELSE
1325 492 : IF (psolver == pw_poisson_implicit) THEN
1326 0 : implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1327 0 : bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1328 0 : SELECT CASE (bc)
1329 : CASE (MIXED_PERIODIC_BC, MIXED_BC)
1330 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1331 0 : "Core Hamiltonian energy: ", energy%core, &
1332 0 : "Hartree energy: ", implicit_ps_ehartree, &
1333 0 : "Electric enthalpy: ", energy%hartree, &
1334 0 : "Exchange-correlation energy: ", energy%exc
1335 : CASE (PERIODIC_BC, NEUMANN_BC)
1336 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1337 0 : "Core Hamiltonian energy: ", energy%core, &
1338 0 : "Hartree energy: ", energy%hartree, &
1339 0 : "Exchange-correlation energy: ", energy%exc
1340 : END SELECT
1341 : ELSE
1342 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1343 492 : "Core Hamiltonian energy: ", energy%core, &
1344 492 : "Hartree energy: ", energy%hartree, &
1345 984 : "Exchange-correlation energy: ", energy%exc
1346 : END IF
1347 : END IF
1348 : END IF
1349 :
1350 492 : IF (dft_control%apply_external_density) THEN
1351 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1352 0 : "Integral of the (density * v_xc): ", energy%exc
1353 : END IF
1354 :
1355 492 : IF (energy%e_hartree /= 0.0_dp) THEN
1356 : WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
1357 460 : "Coulomb (electron-electron) energy: ", energy%e_hartree
1358 : END IF
1359 492 : IF (energy%dispersion /= 0.0_dp) THEN
1360 : WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
1361 0 : "Dispersion energy: ", energy%dispersion
1362 : END IF
1363 492 : IF (energy%efield /= 0.0_dp) THEN
1364 : WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
1365 0 : "Electric field interaction energy: ", energy%efield
1366 : END IF
1367 492 : IF (energy%gcp /= 0.0_dp) THEN
1368 : WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
1369 0 : "gCP energy: ", energy%gcp
1370 : END IF
1371 492 : IF (dft_control%qs_control%gapw) THEN
1372 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1373 32 : "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit, &
1374 64 : "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
1375 : END IF
1376 492 : IF (dft_control%qs_control%gapw_xc) THEN
1377 : WRITE (UNIT=output_unit, FMT="(/,(T3,A,T61,F20.10))") &
1378 0 : "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit
1379 : END IF
1380 492 : IF (dft_control%dft_plus_u) THEN
1381 : WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
1382 0 : "DFT+U energy:", energy%dft_plus_u
1383 : END IF
1384 492 : IF (qs_env%qmmm) THEN
1385 : WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
1386 0 : "QM/MM Electrostatic energy: ", energy%qmmm_el
1387 0 : IF (qs_env%qmmm_env_qm%image_charge) THEN
1388 : WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
1389 0 : "QM/MM image charge energy: ", energy%image_charge
1390 : END IF
1391 : END IF
1392 492 : IF (dft_control%qs_control%mulliken_restraint) THEN
1393 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") &
1394 0 : "Mulliken restraint (order_p,energy) : ", mulliken_order_p, energy%mulliken
1395 : END IF
1396 492 : IF (dft_control%qs_control%ddapc_restraint) THEN
1397 40 : DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
1398 : ddapc_order_p = &
1399 20 : dft_control%qs_control%ddapc_restraint_control(n)%ddapc_order_p
1400 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") &
1401 40 : "DDAPC restraint (order_p,energy) : ", ddapc_order_p, energy%ddapc_restraint(n)
1402 : END DO
1403 : END IF
1404 492 : IF (dft_control%qs_control%s2_restraint) THEN
1405 0 : s2_order_p = dft_control%qs_control%s2_restraint_control%s2_order_p
1406 : WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") &
1407 0 : "S2 restraint (order_p,energy) : ", s2_order_p, energy%s2_restraint
1408 : END IF
1409 492 : IF (energy%core_cneo /= 0.0_dp) THEN
1410 : WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
1411 0 : "CNEO| quantum nuclear core energy: ", energy%core_cneo
1412 : END IF
1413 :
1414 : END IF ! output_unit
1415 : CALL cp_print_key_finished_output(output_unit, logger, input, &
1416 136105 : "DFT%SCF%PRINT%DETAILED_ENERGY")
1417 :
1418 136105 : END SUBROUTINE print_detailed_energy
1419 :
1420 : ! **************************************************************************************************
1421 : !> \brief compute matrix_vxc, defined via the potential created by qs_vxc_create
1422 : !> ignores things like tau functional, gapw, sic, ...
1423 : !> so only OK for GGA & GPW right now
1424 : !> \param qs_env ...
1425 : !> \param v_rspace ...
1426 : !> \param matrix_vxc ...
1427 : !> \param gapw_full_basis ...
1428 : !> \par History
1429 : !> created 23.10.2012 [Joost VandeVondele]
1430 : !> \author
1431 : ! **************************************************************************************************
1432 10 : SUBROUTINE compute_matrix_vxc(qs_env, v_rspace, matrix_vxc, gapw_full_basis)
1433 : TYPE(qs_environment_type), POINTER :: qs_env
1434 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
1435 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
1436 : LOGICAL, INTENT(IN), OPTIONAL :: gapw_full_basis
1437 :
1438 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_matrix_vxc'
1439 :
1440 : INTEGER :: handle, ispin
1441 : LOGICAL :: gapw
1442 10 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1443 : TYPE(dft_control_type), POINTER :: dft_control
1444 :
1445 10 : CALL timeset(routineN, handle)
1446 :
1447 : ! create the matrix using matrix_ks as a template
1448 10 : IF (ASSOCIATED(matrix_vxc)) THEN
1449 0 : CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1450 : END IF
1451 10 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
1452 44 : ALLOCATE (matrix_vxc(SIZE(matrix_ks)))
1453 24 : DO ispin = 1, SIZE(matrix_ks)
1454 14 : NULLIFY (matrix_vxc(ispin)%matrix)
1455 14 : CALL dbcsr_init_p(matrix_vxc(ispin)%matrix)
1456 : CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, &
1457 14 : name="Matrix VXC of spin "//cp_to_string(ispin))
1458 24 : CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp)
1459 : END DO
1460 :
1461 : ! and integrate
1462 10 : CALL get_qs_env(qs_env, dft_control=dft_control)
1463 10 : gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1464 10 : IF (PRESENT(gapw_full_basis)) THEN
1465 0 : IF (gapw_full_basis) gapw = .FALSE.
1466 : END IF
1467 24 : DO ispin = 1, SIZE(matrix_ks)
1468 : CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1469 : hmat=matrix_vxc(ispin), &
1470 : qs_env=qs_env, &
1471 : calculate_forces=.FALSE., &
1472 14 : gapw=gapw)
1473 : ! scale by the volume element... should really become part of integrate_v_rspace
1474 24 : CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw_grid%dvol)
1475 : END DO
1476 :
1477 10 : CALL timestop(handle)
1478 :
1479 10 : END SUBROUTINE compute_matrix_vxc
1480 :
1481 : ! **************************************************************************************************
1482 : !> \brief Build the XC potential matrix for k-point/image-resolved KS matrices.
1483 : !> \param qs_env Quickstep environment
1484 : !> \param v_rspace XC potential on the real-space grid
1485 : !> \param matrix_vxc_kp k-point/image-resolved XC potential matrix
1486 : !> \param gapw_full_basis ...
1487 : !> \author
1488 : ! **************************************************************************************************
1489 0 : SUBROUTINE compute_matrix_vxc_kp(qs_env, v_rspace, matrix_vxc_kp, gapw_full_basis)
1490 : TYPE(qs_environment_type), POINTER :: qs_env
1491 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
1492 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_vxc_kp
1493 : LOGICAL, INTENT(IN), OPTIONAL :: gapw_full_basis
1494 :
1495 : CHARACTER(LEN=*), PARAMETER :: routineN = 'compute_matrix_vxc_kp'
1496 :
1497 : INTEGER :: handle, img, ispin
1498 : LOGICAL :: gapw
1499 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat
1500 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_kp
1501 : TYPE(dft_control_type), POINTER :: dft_control
1502 :
1503 0 : CALL timeset(routineN, handle)
1504 :
1505 : ! create the matrix using matrix_ks as a template
1506 0 : IF (ASSOCIATED(matrix_vxc_kp)) THEN
1507 0 : CALL dbcsr_deallocate_matrix_set(matrix_vxc_kp)
1508 : END IF
1509 0 : CALL get_qs_env(qs_env, dft_control=dft_control, matrix_ks_kp=matrix_ks_kp)
1510 0 : ALLOCATE (matrix_vxc_kp(SIZE(matrix_ks_kp, 1), SIZE(matrix_ks_kp, 2)))
1511 0 : DO img = 1, SIZE(matrix_ks_kp, 2)
1512 0 : DO ispin = 1, SIZE(matrix_ks_kp, 1)
1513 0 : NULLIFY (matrix_vxc_kp(ispin, img)%matrix)
1514 0 : CALL dbcsr_init_p(matrix_vxc_kp(ispin, img)%matrix)
1515 : CALL dbcsr_copy(matrix_vxc_kp(ispin, img)%matrix, matrix_ks_kp(ispin, img)%matrix, &
1516 0 : name="Matrix VXC of spin "//cp_to_string(ispin)//" image "//cp_to_string(img))
1517 0 : CALL dbcsr_set(matrix_vxc_kp(ispin, img)%matrix, 0.0_dp)
1518 : END DO
1519 : END DO
1520 :
1521 : ! and integrate
1522 0 : gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1523 0 : IF (PRESENT(gapw_full_basis)) THEN
1524 0 : IF (gapw_full_basis) gapw = .FALSE.
1525 : END IF
1526 0 : DO ispin = 1, SIZE(matrix_ks_kp, 1)
1527 0 : ksmat => matrix_vxc_kp(ispin, :)
1528 : CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1529 : hmat_kp=ksmat, &
1530 : qs_env=qs_env, &
1531 : calculate_forces=.FALSE., &
1532 0 : gapw=gapw)
1533 : ! scale by the volume element... should really become part of integrate_v_rspace
1534 0 : DO img = 1, SIZE(matrix_ks_kp, 2)
1535 0 : CALL dbcsr_scale(matrix_vxc_kp(ispin, img)%matrix, v_rspace(ispin)%pw_grid%dvol)
1536 : END DO
1537 : END DO
1538 :
1539 0 : CALL timestop(handle)
1540 :
1541 0 : END SUBROUTINE compute_matrix_vxc_kp
1542 :
1543 : ! **************************************************************************************************
1544 : !> \brief Sum up all potentials defined on the grid and integrate
1545 : !>
1546 : !> \param qs_env ...
1547 : !> \param ks_matrix ...
1548 : !> \param rho ...
1549 : !> \param my_rho ...
1550 : !> \param vppl_rspace ...
1551 : !> \param v_rspace_new ...
1552 : !> \param v_rspace_new_aux_fit ...
1553 : !> \param v_tau_rspace ...
1554 : !> \param v_tau_rspace_aux_fit ...
1555 : !> \param v_sic_rspace ...
1556 : !> \param v_spin_ddapc_rest_r ...
1557 : !> \param v_sccs_rspace ...
1558 : !> \param v_rspace_embed ...
1559 : !> \param cdft_control ...
1560 : !> \param calculate_forces ...
1561 : !> \par History
1562 : !> - refactoring 04.03.2011 [MI]
1563 : !> - SCCS implementation (16.10.2013,MK)
1564 : !> \author
1565 : ! **************************************************************************************************
1566 125433 : SUBROUTINE sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, &
1567 : vppl_rspace, v_rspace_new, &
1568 : v_rspace_new_aux_fit, v_tau_rspace, &
1569 : v_tau_rspace_aux_fit, &
1570 : v_sic_rspace, v_spin_ddapc_rest_r, &
1571 : v_sccs_rspace, v_rspace_embed, cdft_control, &
1572 : calculate_forces)
1573 :
1574 : TYPE(qs_environment_type), POINTER :: qs_env
1575 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
1576 : TYPE(qs_rho_type), POINTER :: rho
1577 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: my_rho
1578 : TYPE(pw_r3d_rs_type), POINTER :: vppl_rspace
1579 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_rspace_new_aux_fit, &
1580 : v_tau_rspace, v_tau_rspace_aux_fit
1581 : TYPE(pw_r3d_rs_type), POINTER :: v_sic_rspace, v_spin_ddapc_rest_r, &
1582 : v_sccs_rspace
1583 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_embed
1584 : TYPE(cdft_control_type), POINTER :: cdft_control
1585 : LOGICAL, INTENT(in) :: calculate_forces
1586 :
1587 : CHARACTER(LEN=*), PARAMETER :: routineN = 'sum_up_and_integrate'
1588 :
1589 : CHARACTER(LEN=default_string_length) :: basis_type
1590 : INTEGER :: handle, igroup, ikind, img, ispin, &
1591 : nkind, nspins
1592 125433 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1593 : LOGICAL :: do_ppl, gapw, gapw_composite_direct_ao, &
1594 : gapw_composite_reference, gapw_xc, &
1595 : lrigpw, rigpw, use_work_v_rspace
1596 : REAL(KIND=dp) :: csign, dvol, fadm
1597 : TYPE(admm_type), POINTER :: admm_env
1598 125433 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1599 125433 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao, rho_ao_nokp, smat
1600 125433 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, &
1601 125433 : matrix_ks_aux_fit_dft, rho_ao_aux, &
1602 125433 : rho_ao_kp
1603 : TYPE(dft_control_type), POINTER :: dft_control
1604 : TYPE(kpoint_type), POINTER :: kpoints
1605 : TYPE(lri_density_type), POINTER :: lri_density
1606 : TYPE(lri_environment_type), POINTER :: lri_env
1607 125433 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
1608 : TYPE(mp_para_env_type), POINTER :: para_env
1609 : TYPE(pw_env_type), POINTER :: pw_env
1610 : TYPE(pw_poisson_type), POINTER :: poisson_env
1611 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1612 : TYPE(pw_r3d_rs_type), POINTER :: v_rspace, v_rspace_used, vee
1613 : TYPE(pw_r3d_rs_type), TARGET :: v_rspace_work
1614 : TYPE(qs_ks_env_type), POINTER :: ks_env
1615 : TYPE(qs_rho_type), POINTER :: rho_aux_fit
1616 : TYPE(section_vals_type), POINTER :: input, xc_section
1617 : TYPE(task_list_type), POINTER :: task_list
1618 :
1619 125433 : CALL timeset(routineN, handle)
1620 :
1621 125433 : NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, &
1622 125433 : v_rspace, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, &
1623 125433 : ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, &
1624 125433 : rho_ao_nokp, ks_env, admm_env, task_list, v_rspace_used, input, xc_section)
1625 :
1626 : CALL get_qs_env(qs_env, &
1627 : dft_control=dft_control, &
1628 : input=input, &
1629 : pw_env=pw_env, &
1630 : v_hartree_rspace=v_rspace, &
1631 125433 : vee=vee)
1632 :
1633 125433 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1634 125433 : CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
1635 125433 : gapw = dft_control%qs_control%gapw
1636 125433 : gapw_xc = dft_control%qs_control%gapw_xc
1637 125433 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
1638 : gapw_composite_reference = native_skala_gapw_composite_reference(xc_section) .AND. &
1639 125433 : (gapw .OR. gapw_xc)
1640 : gapw_composite_direct_ao = gapw_composite_reference .AND. &
1641 0 : native_skala_gapw_composite_direct_ao(xc_section)
1642 125433 : do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
1643 :
1644 125433 : rigpw = dft_control%qs_control%rigpw
1645 125433 : lrigpw = dft_control%qs_control%lrigpw
1646 125433 : IF (lrigpw .OR. rigpw) THEN
1647 : CALL get_qs_env(qs_env, &
1648 : lri_env=lri_env, &
1649 : lri_density=lri_density, &
1650 492 : atomic_kind_set=atomic_kind_set)
1651 : END IF
1652 :
1653 125433 : nspins = dft_control%nspins
1654 :
1655 : ! sum up potentials and integrate
1656 125433 : IF (ASSOCIATED(v_rspace_new)) THEN
1657 251541 : DO ispin = 1, nspins
1658 137118 : IF (gapw_composite_reference) THEN
1659 : ! The direct AO diagnostic uses the full ORB basis. The reconstructed path uses
1660 : ! the soft basis here and adds the hard-minus-soft adjoint through update_ks_atom.
1661 0 : CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
1662 0 : rho_ao => rho_ao_kp(ispin, :)
1663 0 : ksmat => ks_matrix(ispin, :)
1664 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1665 : pmat_kp=rho_ao, hmat_kp=ksmat, &
1666 : qs_env=qs_env, &
1667 : calculate_forces=calculate_forces, &
1668 0 : gapw=gapw .AND. .NOT. gapw_composite_direct_ao)
1669 0 : CALL pw_copy(v_rspace, v_rspace_new(ispin))
1670 137118 : ELSE IF (gapw_xc) THEN
1671 : ! SIC not implemented (or at least not tested)
1672 4450 : CPASSERT(dft_control%sic_method_id == sic_none)
1673 : !Only the xc potential, because it has to be integrated with the soft basis
1674 4450 : CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
1675 :
1676 : ! add the xc part due to v_rspace soft
1677 4450 : rho_ao => rho_ao_kp(ispin, :)
1678 4450 : ksmat => ks_matrix(ispin, :)
1679 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1680 : pmat_kp=rho_ao, hmat_kp=ksmat, &
1681 : qs_env=qs_env, &
1682 : calculate_forces=calculate_forces, &
1683 4450 : gapw=gapw_xc)
1684 :
1685 : ! Now the Hartree potential to be integrated with the full basis
1686 4450 : CALL pw_copy(v_rspace, v_rspace_new(ispin))
1687 : ELSE
1688 : ! Add v_hartree + v_xc = v_rspace_new
1689 132668 : CALL pw_axpy(v_rspace, v_rspace_new(ispin), 1.0_dp, v_rspace_new(ispin)%pw_grid%dvol)
1690 : END IF ! gapw_xc
1691 137118 : IF (dft_control%qs_control%ddapc_explicit_potential) THEN
1692 184 : IF (dft_control%qs_control%ddapc_restraint_is_spin) THEN
1693 184 : IF (ispin == 1) THEN
1694 92 : CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1695 : ELSE
1696 92 : CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), -1.0_dp)
1697 : END IF
1698 : ELSE
1699 0 : CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1700 : END IF
1701 : END IF
1702 : ! CDFT constraint contribution
1703 137118 : IF (dft_control%qs_control%cdft) THEN
1704 12392 : DO igroup = 1, SIZE(cdft_control%group)
1705 7216 : SELECT CASE (cdft_control%group(igroup)%constraint_type)
1706 : CASE (cdft_charge_constraint)
1707 32 : csign = 1.0_dp
1708 : CASE (cdft_magnetization_constraint)
1709 32 : IF (ispin == 1) THEN
1710 : csign = 1.0_dp
1711 : ELSE
1712 16 : csign = -1.0_dp
1713 : END IF
1714 : CASE (cdft_alpha_constraint)
1715 2016 : csign = 1.0_dp
1716 2016 : IF (ispin == 2) CYCLE
1717 : CASE (cdft_beta_constraint)
1718 2016 : csign = 1.0_dp
1719 2016 : IF (ispin == 1) CYCLE
1720 : CASE DEFAULT
1721 7216 : CPABORT("Unknown constraint type.")
1722 : END SELECT
1723 : CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_new(ispin), &
1724 12392 : csign*cdft_control%strength(igroup))
1725 : END DO
1726 : END IF
1727 : ! functional derivative of the Hartree energy wrt the density in the presence of dielectric
1728 : ! (vhartree + v_eps); v_eps is nonzero only if the dielectric constant is defind as a function
1729 : ! of the charge density
1730 137118 : IF (poisson_env%parameters%solver == pw_poisson_implicit) THEN
1731 486 : dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1732 486 : CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_new(ispin), dvol)
1733 : END IF
1734 : ! Add SCCS contribution
1735 137118 : IF (dft_control%do_sccs) THEN
1736 282 : CALL pw_axpy(v_sccs_rspace, v_rspace_new(ispin))
1737 : END IF
1738 : ! External electrostatic potential
1739 137118 : IF (dft_control%apply_external_potential) THEN
1740 : CALL qmmm_modify_hartree_pot(v_hartree=v_rspace_new(ispin), &
1741 364 : v_qmmm=vee, scale=-1.0_dp)
1742 : END IF
1743 137118 : IF (do_ppl) THEN
1744 66 : CPASSERT(.NOT. gapw)
1745 66 : CALL pw_axpy(vppl_rspace, v_rspace_new(ispin), vppl_rspace%pw_grid%dvol)
1746 : END IF
1747 : ! the electrostatic sic contribution
1748 137238 : SELECT CASE (dft_control%sic_method_id)
1749 : CASE (sic_none)
1750 : !
1751 : CASE (sic_mauri_us, sic_mauri_spz)
1752 120 : IF (ispin == 1) THEN
1753 60 : CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1754 : ELSE
1755 60 : CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), 1.0_dp)
1756 : END IF
1757 : CASE (sic_ad)
1758 137118 : CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1759 : CASE (sic_eo)
1760 : ! NOTHING TO BE DONE
1761 : END SELECT
1762 : ! DFT embedding
1763 137118 : IF (dft_control%apply_embed_pot) THEN
1764 838 : CALL pw_axpy(v_rspace_embed(ispin), v_rspace_new(ispin), v_rspace_embed(ispin)%pw_grid%dvol)
1765 838 : CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1766 : END IF
1767 137118 : IF (lrigpw) THEN
1768 478 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1769 478 : CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1770 1430 : DO ikind = 1, nkind
1771 308200 : lri_v_int(ikind)%v_int = 0.0_dp
1772 : END DO
1773 : CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1774 478 : lri_v_int, calculate_forces, "LRI_AUX")
1775 1430 : DO ikind = 1, nkind
1776 614970 : CALL para_env%sum(lri_v_int(ikind)%v_int)
1777 : END DO
1778 478 : IF (lri_env%exact_1c_terms) THEN
1779 36 : rho_ao => my_rho(ispin, :)
1780 36 : ksmat => ks_matrix(ispin, :)
1781 : CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), ksmat(1)%matrix, &
1782 : rho_ao(1)%matrix, qs_env, &
1783 36 : calculate_forces, "ORB")
1784 : END IF
1785 478 : IF (lri_env%ppl_ri) THEN
1786 8 : CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1787 : END IF
1788 136640 : ELSE IF (rigpw) THEN
1789 26 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1790 26 : CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1791 52 : DO ikind = 1, nkind
1792 1144 : lri_v_int(ikind)%v_int = 0.0_dp
1793 : END DO
1794 : CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1795 26 : lri_v_int, calculate_forces, "RI_HXC")
1796 52 : DO ikind = 1, nkind
1797 2236 : CALL para_env%sum(lri_v_int(ikind)%v_int)
1798 : END DO
1799 : ELSE
1800 136614 : rho_ao => my_rho(ispin, :)
1801 136614 : ksmat => ks_matrix(ispin, :)
1802 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1803 : pmat_kp=rho_ao, hmat_kp=ksmat, &
1804 : qs_env=qs_env, &
1805 : calculate_forces=calculate_forces, &
1806 136614 : gapw=gapw .AND. .NOT. gapw_composite_direct_ao)
1807 : END IF
1808 251541 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
1809 : END DO ! ispin
1810 :
1811 114483 : SELECT CASE (dft_control%sic_method_id)
1812 : CASE (sic_none)
1813 : CASE (sic_mauri_us, sic_mauri_spz, sic_ad)
1814 60 : CALL auxbas_pw_pool%give_back_pw(v_sic_rspace)
1815 114483 : DEALLOCATE (v_sic_rspace)
1816 : END SELECT
1817 114423 : DEALLOCATE (v_rspace_new)
1818 :
1819 : ELSE
1820 : ! not implemented (or at least not tested)
1821 11010 : CPASSERT(dft_control%sic_method_id == sic_none)
1822 11010 : CPASSERT(.NOT. dft_control%qs_control%ddapc_restraint_is_spin)
1823 24306 : DO ispin = 1, nspins
1824 13296 : use_work_v_rspace = dft_control%qs_control%cdft
1825 13296 : IF (use_work_v_rspace) THEN
1826 168 : CALL auxbas_pw_pool%create_pw(v_rspace_work)
1827 168 : CALL pw_copy(v_rspace, v_rspace_work)
1828 168 : v_rspace_used => v_rspace_work
1829 : ELSE
1830 13128 : v_rspace_used => v_rspace
1831 : END IF
1832 : ! CDFT constraint contribution
1833 13296 : IF (dft_control%qs_control%cdft) THEN
1834 336 : DO igroup = 1, SIZE(cdft_control%group)
1835 168 : SELECT CASE (cdft_control%group(igroup)%constraint_type)
1836 : CASE (cdft_charge_constraint)
1837 0 : csign = 1.0_dp
1838 : CASE (cdft_magnetization_constraint)
1839 0 : IF (ispin == 1) THEN
1840 : csign = 1.0_dp
1841 : ELSE
1842 0 : csign = -1.0_dp
1843 : END IF
1844 : CASE (cdft_alpha_constraint)
1845 0 : csign = 1.0_dp
1846 0 : IF (ispin == 2) CYCLE
1847 : CASE (cdft_beta_constraint)
1848 0 : csign = 1.0_dp
1849 0 : IF (ispin == 1) CYCLE
1850 : CASE DEFAULT
1851 168 : CPABORT("Unknown constraint type.")
1852 : END SELECT
1853 : CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_used, &
1854 336 : csign*cdft_control%strength(igroup))
1855 : END DO
1856 : END IF
1857 : ! extra contribution attributed to the dependency of the dielectric constant to the charge density
1858 13296 : IF (poisson_env%parameters%solver == pw_poisson_implicit) THEN
1859 0 : dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1860 0 : CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_used, dvol)
1861 : END IF
1862 : ! Add SCCS contribution
1863 13296 : IF (dft_control%do_sccs) THEN
1864 0 : CALL pw_axpy(v_sccs_rspace, v_rspace_used)
1865 : END IF
1866 : ! DFT embedding
1867 13296 : IF (dft_control%apply_embed_pot) THEN
1868 12 : CALL pw_axpy(v_rspace_embed(ispin), v_rspace_used, v_rspace_embed(ispin)%pw_grid%dvol)
1869 12 : CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1870 : END IF
1871 13296 : IF (lrigpw) THEN
1872 0 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1873 0 : CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1874 0 : DO ikind = 1, nkind
1875 0 : lri_v_int(ikind)%v_int = 0.0_dp
1876 : END DO
1877 : CALL integrate_v_rspace_one_center(v_rspace_used, qs_env, &
1878 0 : lri_v_int, calculate_forces, "LRI_AUX")
1879 0 : DO ikind = 1, nkind
1880 0 : CALL para_env%sum(lri_v_int(ikind)%v_int)
1881 : END DO
1882 0 : IF (lri_env%exact_1c_terms) THEN
1883 0 : rho_ao => my_rho(ispin, :)
1884 0 : ksmat => ks_matrix(ispin, :)
1885 : CALL integrate_v_rspace_diagonal(v_rspace_used, ksmat(1)%matrix, &
1886 : rho_ao(1)%matrix, qs_env, &
1887 0 : calculate_forces, "ORB")
1888 : END IF
1889 0 : IF (lri_env%ppl_ri) THEN
1890 0 : CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1891 : END IF
1892 13296 : ELSE IF (rigpw) THEN
1893 0 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1894 0 : CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1895 0 : DO ikind = 1, nkind
1896 0 : lri_v_int(ikind)%v_int = 0.0_dp
1897 : END DO
1898 : CALL integrate_v_rspace_one_center(v_rspace_used, qs_env, &
1899 0 : lri_v_int, calculate_forces, "RI_HXC")
1900 0 : DO ikind = 1, nkind
1901 0 : CALL para_env%sum(lri_v_int(ikind)%v_int)
1902 : END DO
1903 : ELSE
1904 13296 : rho_ao => my_rho(ispin, :)
1905 13296 : ksmat => ks_matrix(ispin, :)
1906 : CALL integrate_v_rspace(v_rspace=v_rspace_used, &
1907 : pmat_kp=rho_ao, &
1908 : hmat_kp=ksmat, &
1909 : qs_env=qs_env, &
1910 : calculate_forces=calculate_forces, &
1911 13296 : gapw=gapw)
1912 : END IF
1913 24306 : IF (use_work_v_rspace) CALL auxbas_pw_pool%give_back_pw(v_rspace_work)
1914 : END DO
1915 : END IF ! ASSOCIATED(v_rspace_new)
1916 :
1917 : ! **** LRIGPW: KS matrix from integrated potential
1918 125433 : IF (lrigpw) THEN
1919 466 : CALL get_qs_env(qs_env, ks_env=ks_env)
1920 466 : CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1921 466 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1922 944 : DO ispin = 1, nspins
1923 478 : ksmat => ks_matrix(ispin, :)
1924 : CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set, &
1925 944 : cell_to_index=cell_to_index)
1926 : END DO
1927 466 : IF (calculate_forces) THEN
1928 24 : CALL calculate_lri_forces(lri_env, lri_density, qs_env, rho_ao_kp, atomic_kind_set)
1929 : END IF
1930 124967 : ELSE IF (rigpw) THEN
1931 26 : CALL get_qs_env(qs_env, matrix_s=smat)
1932 52 : DO ispin = 1, nspins
1933 : CALL calculate_ri_ks_matrix(lri_env, lri_v_int, ks_matrix(ispin, 1)%matrix, &
1934 52 : smat(1)%matrix, atomic_kind_set, ispin)
1935 : END DO
1936 26 : IF (calculate_forces) THEN
1937 2 : rho_ao_nokp => rho_ao_kp(:, 1)
1938 2 : CALL calculate_ri_forces(lri_env, lri_density, qs_env, rho_ao_nokp, atomic_kind_set)
1939 : END IF
1940 : END IF
1941 :
1942 125433 : IF (ASSOCIATED(v_tau_rspace)) THEN
1943 3332 : IF (lrigpw .OR. rigpw) THEN
1944 0 : CPABORT("LRIGPW/RIGPW not implemented for meta-GGAs")
1945 : END IF
1946 7260 : DO ispin = 1, nspins
1947 3928 : CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1948 :
1949 3928 : rho_ao => rho_ao_kp(ispin, :)
1950 3928 : ksmat => ks_matrix(ispin, :)
1951 : CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1952 : pmat_kp=rho_ao, hmat_kp=ksmat, &
1953 : qs_env=qs_env, &
1954 : calculate_forces=calculate_forces, compute_tau=.TRUE., &
1955 : gapw=(gapw .OR. gapw_xc) .AND. &
1956 6956 : .NOT. gapw_composite_direct_ao)
1957 7260 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1958 : END DO
1959 3332 : DEALLOCATE (v_tau_rspace)
1960 : END IF
1961 :
1962 : ! Add contributions from ADMM if requested
1963 125433 : IF (dft_control%do_admm) THEN
1964 12318 : CALL get_qs_env(qs_env, admm_env=admm_env)
1965 : CALL get_admm_env(admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
1966 12318 : matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft)
1967 12318 : CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
1968 12318 : IF (ASSOCIATED(v_rspace_new_aux_fit)) THEN
1969 18806 : DO ispin = 1, nspins
1970 : ! Calculate the xc potential
1971 10262 : CALL pw_scale(v_rspace_new_aux_fit(ispin), v_rspace_new_aux_fit(ispin)%pw_grid%dvol)
1972 :
1973 : ! set matrix_ks_aux_fit_dft = matrix_ks_aux_fit(k_HF)
1974 25464 : DO img = 1, dft_control%nimages
1975 : CALL dbcsr_copy(matrix_ks_aux_fit_dft(ispin, img)%matrix, matrix_ks_aux_fit(ispin, img)%matrix, &
1976 25464 : name="DFT exch. part of matrix_ks_aux_fit")
1977 : END DO
1978 :
1979 : ! Add potential to ks_matrix aux_fit, skip integration if no DFT correction
1980 :
1981 10262 : IF (admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
1982 :
1983 : !GPW by default. IF GAPW, then take relevant task list and basis
1984 10262 : CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
1985 10262 : basis_type = "AUX_FIT"
1986 10262 : IF (admm_env%do_gapw) THEN
1987 3832 : task_list => admm_env%admm_gapw_env%task_list
1988 3832 : basis_type = "AUX_FIT_SOFT"
1989 : END IF
1990 10262 : fadm = 1.0_dp
1991 : ! Calculate bare scaling of force according to Merlot, 1. IF: ADMMP, 2. IF: ADMMS,
1992 10262 : IF (admm_env%do_admmp) THEN
1993 444 : fadm = admm_env%gsi(ispin)**2
1994 9818 : ELSE IF (admm_env%do_admms) THEN
1995 476 : fadm = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)
1996 : END IF
1997 :
1998 10262 : rho_ao => rho_ao_aux(ispin, :)
1999 10262 : ksmat => matrix_ks_aux_fit(ispin, :)
2000 :
2001 : CALL integrate_v_rspace(v_rspace=v_rspace_new_aux_fit(ispin), &
2002 : pmat_kp=rho_ao, &
2003 : hmat_kp=ksmat, &
2004 : qs_env=qs_env, &
2005 : calculate_forces=calculate_forces, &
2006 : force_adm=fadm, &
2007 : gapw=.FALSE., & !even if actual GAPW calculation, want to use AUX_FIT_SOFT
2008 : basis_type=basis_type, &
2009 10262 : task_list_external=task_list)
2010 : END IF
2011 :
2012 : ! matrix_ks_aux_fit_dft(x_DFT)=matrix_ks_aux_fit_dft(old,k_HF)-matrix_ks_aux_fit(k_HF-x_DFT)
2013 25464 : DO img = 1, dft_control%nimages
2014 : CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, &
2015 25464 : matrix_ks_aux_fit(ispin, img)%matrix, 1.0_dp, -1.0_dp)
2016 : END DO
2017 :
2018 18806 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new_aux_fit(ispin))
2019 : END DO
2020 8544 : DEALLOCATE (v_rspace_new_aux_fit)
2021 : END IF
2022 : ! Clean up v_tau_rspace_aux_fit, which is actually not needed
2023 12318 : IF (ASSOCIATED(v_tau_rspace_aux_fit)) THEN
2024 0 : DO ispin = 1, nspins
2025 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace_aux_fit(ispin))
2026 : END DO
2027 0 : DEALLOCATE (v_tau_rspace_aux_fit)
2028 : END IF
2029 : END IF
2030 :
2031 125433 : IF (dft_control%apply_embed_pot) DEALLOCATE (v_rspace_embed)
2032 :
2033 125433 : CALL timestop(handle)
2034 :
2035 125433 : END SUBROUTINE sum_up_and_integrate
2036 :
2037 : !**************************************************************************
2038 : !> \brief Calculate the ZMP potential and energy as in Zhao, Morrison Parr
2039 : !> PRA 50i, 2138 (1994)
2040 : !> V_c^\lambda defined as int_rho-rho_0/r-r' or rho-rho_0 times a Lagrange
2041 : !> multiplier, plus Fermi-Amaldi potential that should give the V_xc in the
2042 : !> limit \lambda --> \infty
2043 : !>
2044 : !> \param qs_env ...
2045 : !> \param v_rspace_new ...
2046 : !> \param rho ...
2047 : !> \param exc ...
2048 : !> \author D. Varsano [daniele.varsano@nano.cnr.it]
2049 : ! **************************************************************************************************
2050 106 : SUBROUTINE calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
2051 :
2052 : TYPE(qs_environment_type), POINTER :: qs_env
2053 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new
2054 : TYPE(qs_rho_type), POINTER :: rho
2055 : REAL(KIND=dp) :: exc
2056 :
2057 : CHARACTER(*), PARAMETER :: routineN = 'calculate_zmp_potential'
2058 :
2059 : INTEGER :: handle, my_val, nelectron, nspins
2060 : INTEGER, DIMENSION(2) :: nelectron_spin
2061 : LOGICAL :: do_zmp_read, fermi_amaldi
2062 : REAL(KIND=dp) :: lambda
2063 106 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_ext_r
2064 : TYPE(dft_control_type), POINTER :: dft_control
2065 106 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ext_g, rho_g
2066 : TYPE(pw_env_type), POINTER :: pw_env
2067 : TYPE(pw_poisson_type), POINTER :: poisson_env
2068 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2069 : TYPE(pw_r3d_rs_type) :: v_xc_rspace
2070 106 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
2071 : TYPE(qs_ks_env_type), POINTER :: ks_env
2072 : TYPE(section_vals_type), POINTER :: ext_den_section, input
2073 :
2074 : !, v_h_gspace, &
2075 :
2076 106 : CALL timeset(routineN, handle)
2077 106 : NULLIFY (auxbas_pw_pool)
2078 106 : NULLIFY (pw_env)
2079 106 : NULLIFY (poisson_env)
2080 106 : NULLIFY (v_rspace_new)
2081 106 : NULLIFY (dft_control)
2082 106 : NULLIFY (rho_r, rho_g, tot_rho_ext_r, rho_ext_g)
2083 : CALL get_qs_env(qs_env=qs_env, &
2084 : pw_env=pw_env, &
2085 : ks_env=ks_env, &
2086 : rho=rho, &
2087 : input=input, &
2088 : nelectron_spin=nelectron_spin, &
2089 106 : dft_control=dft_control)
2090 : CALL pw_env_get(pw_env=pw_env, &
2091 : auxbas_pw_pool=auxbas_pw_pool, &
2092 106 : poisson_env=poisson_env)
2093 106 : CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
2094 106 : nspins = 1
2095 212 : ALLOCATE (v_rspace_new(nspins))
2096 106 : CALL auxbas_pw_pool%create_pw(pw=v_rspace_new(1))
2097 106 : CALL auxbas_pw_pool%create_pw(pw=v_xc_rspace)
2098 :
2099 106 : CALL pw_zero(v_rspace_new(1))
2100 106 : do_zmp_read = dft_control%apply_external_vxc
2101 106 : IF (do_zmp_read) THEN
2102 0 : CALL pw_copy(qs_env%external_vxc, v_rspace_new(1))
2103 : exc = accurate_dot_product(v_rspace_new(1)%array, rho_r(1)%array)* &
2104 0 : v_rspace_new(1)%pw_grid%dvol
2105 : ELSE
2106 318 : BLOCK
2107 : REAL(KIND=dp) :: factor
2108 : TYPE(pw_c1d_gs_type) :: rho_eff_gspace, v_xc_gspace
2109 106 : CALL auxbas_pw_pool%create_pw(pw=rho_eff_gspace)
2110 106 : CALL auxbas_pw_pool%create_pw(pw=v_xc_gspace)
2111 106 : CALL pw_zero(rho_eff_gspace)
2112 106 : CALL pw_zero(v_xc_gspace)
2113 106 : CALL pw_zero(v_xc_rspace)
2114 106 : factor = pw_integrate_function(rho_g(1))
2115 : CALL qs_rho_get(qs_env%rho_external, &
2116 : rho_g=rho_ext_g, &
2117 106 : tot_rho_r=tot_rho_ext_r)
2118 106 : factor = tot_rho_ext_r(1)/factor
2119 :
2120 106 : CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
2121 106 : CALL pw_axpy(rho_ext_g(1), rho_eff_gspace, alpha=-1.0_dp)
2122 106 : ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
2123 106 : CALL section_vals_val_get(ext_den_section, "LAMBDA", r_val=lambda)
2124 106 : CALL section_vals_val_get(ext_den_section, "ZMP_CONSTRAINT", i_val=my_val)
2125 106 : CALL section_vals_val_get(ext_den_section, "FERMI_AMALDI", l_val=fermi_amaldi)
2126 :
2127 106 : CALL pw_scale(rho_eff_gspace, a=lambda)
2128 106 : nelectron = nelectron_spin(1)
2129 106 : factor = -1.0_dp/nelectron
2130 106 : CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
2131 :
2132 106 : CALL pw_poisson_solve(poisson_env, rho_eff_gspace, vhartree=v_xc_gspace)
2133 106 : CALL pw_transfer(v_xc_gspace, v_rspace_new(1))
2134 106 : CALL pw_copy(v_rspace_new(1), v_xc_rspace)
2135 :
2136 106 : exc = 0.0_dp
2137 106 : exc = pw_integral_ab(v_rspace_new(1), rho_r(1))
2138 :
2139 : !Note that this is not the xc energy but \int(\rho*v_xc)
2140 : !Vxc---> v_rspace_new
2141 : !Exc---> energy%exc
2142 106 : CALL auxbas_pw_pool%give_back_pw(rho_eff_gspace)
2143 212 : CALL auxbas_pw_pool%give_back_pw(v_xc_gspace)
2144 : END BLOCK
2145 : END IF
2146 :
2147 106 : CALL auxbas_pw_pool%give_back_pw(v_xc_rspace)
2148 :
2149 106 : CALL timestop(handle)
2150 :
2151 106 : END SUBROUTINE calculate_zmp_potential
2152 :
2153 : ! **************************************************************************************************
2154 : !> \brief ...
2155 : !> \param qs_env ...
2156 : !> \param rho ...
2157 : !> \param v_rspace_embed ...
2158 : !> \param dft_control ...
2159 : !> \param embed_corr ...
2160 : !> \param just_energy ...
2161 : ! **************************************************************************************************
2162 822 : SUBROUTINE get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, &
2163 : just_energy)
2164 : TYPE(qs_environment_type), POINTER :: qs_env
2165 : TYPE(qs_rho_type), POINTER :: rho
2166 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_embed
2167 : TYPE(dft_control_type), POINTER :: dft_control
2168 : REAL(KIND=dp) :: embed_corr
2169 : LOGICAL :: just_energy
2170 :
2171 : CHARACTER(*), PARAMETER :: routineN = 'get_embed_potential_energy'
2172 :
2173 : INTEGER :: handle, ispin
2174 : REAL(KIND=dp) :: embed_corr_local
2175 : TYPE(pw_env_type), POINTER :: pw_env
2176 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
2177 822 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
2178 :
2179 822 : CALL timeset(routineN, handle)
2180 :
2181 822 : NULLIFY (auxbas_pw_pool)
2182 822 : NULLIFY (pw_env)
2183 822 : NULLIFY (rho_r)
2184 : CALL get_qs_env(qs_env=qs_env, &
2185 : pw_env=pw_env, &
2186 822 : rho=rho)
2187 : CALL pw_env_get(pw_env=pw_env, &
2188 822 : auxbas_pw_pool=auxbas_pw_pool)
2189 822 : CALL qs_rho_get(rho, rho_r=rho_r)
2190 3722 : ALLOCATE (v_rspace_embed(dft_control%nspins))
2191 :
2192 822 : embed_corr = 0.0_dp
2193 :
2194 2078 : DO ispin = 1, dft_control%nspins
2195 1256 : CALL auxbas_pw_pool%create_pw(pw=v_rspace_embed(ispin))
2196 1256 : CALL pw_zero(v_rspace_embed(ispin))
2197 :
2198 1256 : CALL pw_copy(qs_env%embed_pot, v_rspace_embed(ispin))
2199 1256 : embed_corr_local = 0.0_dp
2200 :
2201 : ! Spin embedding potential in open-shell case
2202 1256 : IF (dft_control%nspins == 2) THEN
2203 868 : IF (ispin == 1) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), 1.0_dp)
2204 868 : IF (ispin == 2) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), -1.0_dp)
2205 : END IF
2206 : ! Integrate the density*potential
2207 1256 : embed_corr_local = pw_integral_ab(v_rspace_embed(ispin), rho_r(ispin))
2208 :
2209 2078 : embed_corr = embed_corr + embed_corr_local
2210 :
2211 : END DO
2212 :
2213 : ! If only energy requiested we delete the potential
2214 822 : IF (just_energy) THEN
2215 692 : DO ispin = 1, dft_control%nspins
2216 692 : CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
2217 : END DO
2218 286 : DEALLOCATE (v_rspace_embed)
2219 : END IF
2220 :
2221 822 : CALL timestop(handle)
2222 :
2223 822 : END SUBROUTINE get_embed_potential_energy
2224 :
2225 : END MODULE qs_ks_utils
|
