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