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