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