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