(git:e7e05ae)
qs_ks_utils.F
Go to the documentation of this file.
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_operations, ONLY: copy_dbcsr_to_fm,&
29  copy_fm_to_dbcsr,&
30  cp_dbcsr_plus_fm_fm_t,&
31  cp_dbcsr_sm_fm_multiply,&
32  dbcsr_allocate_matrix_set,&
33  dbcsr_deallocate_matrix_set
34  USE cp_ddapc, ONLY: cp_ddapc_apply_cd
35  USE cp_fm_struct, ONLY: cp_fm_struct_create,&
36  cp_fm_struct_release,&
37  cp_fm_struct_type
38  USE cp_fm_types, ONLY: cp_fm_create,&
39  cp_fm_get_info,&
40  cp_fm_release,&
41  cp_fm_set_all,&
42  cp_fm_to_fm,&
43  cp_fm_type
44  USE cp_log_handling, ONLY: cp_get_default_logger,&
45  cp_logger_type,&
46  cp_to_string
47  USE cp_output_handling, ONLY: cp_p_file,&
48  cp_print_key_finished_output,&
49  cp_print_key_should_output,&
50  cp_print_key_unit_nr
51  USE dbcsr_api, ONLY: &
52  dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_dot, dbcsr_get_info, dbcsr_init_p, &
53  dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_scale, dbcsr_scale_by_vector, &
54  dbcsr_set, dbcsr_type
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  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  INTEGER, DIMENSION(:), POINTER :: ivec
161  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  REAL(kind=dp), DIMENSION(:), POINTER :: energy_scaling, rvec, scaling
167  TYPE(cell_type), POINTER :: cell
168  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_hfx, matrix_p, mdummy, &
169  mo_derivs, rho_ao
170  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p2
171  TYPE(dbcsr_type), POINTER :: dbcsr_deriv, fm_deriv, fm_scaled, &
172  mo_coeff
173  TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
174  TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
175  TYPE(mp_para_env_type), POINTER :: para_env
176  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  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  IF (.NOT. dft_control%low_spin_roks) RETURN
188 
189  CALL timeset(routinen, handle)
190 
191  NULLIFY (ks_env, rho_ao)
192 
193  ! Test for not compatible options
194  IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
195  CALL cp_abort(__location__, "GAPW/GAPW_XC are not compatible with low spin ROKS method.")
196  END IF
197  IF (dft_control%do_admm) THEN
198  CALL cp_abort(__location__, "ADMM not compatible with low spin ROKS method.")
199  END IF
200  IF (dft_control%do_admm) THEN
201  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  "not compatible with low spin ROKS method.")
204  END IF
205  END IF
206  IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
207  dft_control%qs_control%xtb) THEN
208  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  virial=virial)
220 
221  CALL qs_rho_get(rho, rho_ao=rho_ao)
222 
223  compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
224  xc_section => section_vals_get_subs_vals(input, "DFT%XC")
225  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  cpassert(SIZE(mo_array, 1) == 2)
230  nspin = 2
231  ! we want uniform occupations
232  CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
233  cpassert(uniform_occupation)
234  CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation)
235  cpassert(uniform_occupation)
236  IF (do_hfx .AND. calculate_forces .AND. compute_virial) THEN
237  CALL cp_abort(__location__, "ROKS virial with HFX not available.")
238  END IF
239 
240  NULLIFY (dbcsr_deriv)
241  CALL dbcsr_init_p(dbcsr_deriv)
242  CALL dbcsr_copy(dbcsr_deriv, mo_derivs(1)%matrix)
243  CALL dbcsr_set(dbcsr_deriv, 0.0_dp)
244 
245  ! basic info
246  CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
247  CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_alpha)
248  CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff)
249  CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_beta)
250 
251  ! read the input
252  low_spin_roks_section => section_vals_get_subs_vals(input, "DFT%LOW_SPIN_ROKS")
253 
254  CALL section_vals_val_get(low_spin_roks_section, "ENERGY_SCALING", r_vals=rvec)
255  nterms = SIZE(rvec)
256  ALLOCATE (energy_scaling(nterms))
257  energy_scaling = rvec !? just wondering, should this add up to 1, in which case we should cpp?
258 
259  CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", n_rep_val=n_rep)
260  cpassert(n_rep == nterms)
261  CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=1, i_vals=ivec)
262  nelectron = SIZE(ivec)
263  cpassert(nelectron == k_alpha - k_beta)
264  ALLOCATE (occupations(2, nelectron, nterms))
265  occupations = 0
266  DO iterm = 1, nterms
267  CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=iterm, i_vals=ivec)
268  cpassert(nelectron == SIZE(ivec))
269  in_range = all(ivec >= 1) .AND. all(ivec <= 2)
270  cpassert(in_range)
271  DO k = 1, nelectron
272  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  NULLIFY (matrix_p)
280  CALL dbcsr_allocate_matrix_set(matrix_p, nspin)
281  DO ispin = 1, nspin
282  ALLOCATE (matrix_p(ispin)%matrix)
283  CALL dbcsr_copy(matrix_p(ispin)%matrix, rho_ao(1)%matrix, &
284  name="density matrix low spin roks")
285  CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
286  END DO
287 
288  NULLIFY (matrix_h)
289  CALL dbcsr_allocate_matrix_set(matrix_h, nspin)
290  DO ispin = 1, nspin
291  ALLOCATE (matrix_h(ispin)%matrix)
292  CALL dbcsr_copy(matrix_h(ispin)%matrix, rho_ao(1)%matrix, &
293  name="KS matrix low spin roks")
294  CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
295  END DO
296 
297  IF (do_hfx) THEN
298  NULLIFY (matrix_hfx)
299  CALL dbcsr_allocate_matrix_set(matrix_hfx, nspin)
300  DO ispin = 1, nspin
301  ALLOCATE (matrix_hfx(ispin)%matrix)
302  CALL dbcsr_copy(matrix_hfx(ispin)%matrix, rho_ao(1)%matrix, &
303  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  NULLIFY (tau, vxc_tau, vxc)
310  CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
311 
312  ALLOCATE (rho_r(nspin))
313  ALLOCATE (rho_g(nspin))
314  DO ispin = 1, nspin
315  CALL auxbas_pw_pool%create_pw(rho_r(ispin))
316  CALL auxbas_pw_pool%create_pw(rho_g(ispin))
317  END DO
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  CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
323  NULLIFY (fm_scaled, fm_deriv)
324  CALL dbcsr_init_p(fm_scaled)
325  CALL dbcsr_init_p(fm_deriv)
326  CALL dbcsr_copy(fm_scaled, mo_coeff)
327  CALL dbcsr_copy(fm_deriv, mo_coeff)
328 
329  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  DO iterm = 1, nterms
333 
334  DO ispin = 1, nspin
335  ! compute the proper density matrices with the required occupations
336  CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
337  scaling = 1.0_dp
338  scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
339  CALL dbcsr_copy(fm_scaled, mo_coeff)
340  CALL dbcsr_scale_by_vector(fm_scaled, scaling, side='right')
341  CALL dbcsr_multiply('n', 't', 1.0_dp, mo_coeff, fm_scaled, &
342  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  ks_env=ks_env)
347  END DO
348 
349  ! compute the exchange energies / potential if needed
350  IF (just_energy) THEN
351  exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
352  pw_pool=xc_pw_pool)
353  ELSE
354  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  pw_pool=xc_pw_pool, compute_virial=.false., virial_xc=virial_xc_tmp)
358  END IF
359 
360  energy%exc = energy%exc + energy_scaling(iterm)*exc
361 
362  IF (do_hfx) THEN
363  ! Add Hartree-Fock contribution
364  DO ispin = 1, nspin
365  CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
366  END DO
367  ehfx = energy%ex
368  CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, &
369  recalc_integrals=.false., update_energy=.true.)
370  energy%ex = ehfx + energy_scaling(iterm)*energy%ex
371  END IF
372 
373  ! add the corresponding derivatives to the MO derivatives
374  IF (.NOT. just_energy) THEN
375  ! get the potential in matrix form
376  DO ispin = 1, nspin
377  CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
378  ! use a work_v_rspace
379  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  qs_env=qs_env, calculate_forces=calculate_forces)
382  CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
383  END DO
384  DEALLOCATE (vxc)
385 
386  IF (do_hfx) THEN
387  ! add HFX contribution
388  DO ispin = 1, nspin
389  CALL dbcsr_add(matrix_h(ispin)%matrix, matrix_hfx(ispin)%matrix, &
390  1.0_dp, energy_scaling(iterm))
391  END DO
392  IF (calculate_forces) THEN
393  CALL get_qs_env(qs_env, x_data=x_data, para_env=para_env)
394  IF (x_data(1, 1)%n_rep_hf /= 1) THEN
395  CALL cp_abort(__location__, "Multiple HFX section forces not compatible "// &
396  "with low spin ROKS method.")
397  END IF
398  IF (x_data(1, 1)%do_hfx_ri) THEN
399  CALL cp_abort(__location__, "HFX_RI forces not compatible with low spin ROKS method.")
400  ELSE
401  irep = 1
402  NULLIFY (mdummy)
403  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  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  DO ispin = 1, nspin
414  CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, &
415  0.0_dp, dbcsr_deriv, last_column=k_alpha)
416 
417  scaling = 1.0_dp
418  scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
419  CALL dbcsr_scale_by_vector(dbcsr_deriv, scaling, side='right')
420  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  DO ispin = 1, nspin
429  CALL auxbas_pw_pool%give_back_pw(rho_r(ispin))
430  CALL auxbas_pw_pool%give_back_pw(rho_g(ispin))
431  END DO
432  DEALLOCATE (rho_r, rho_g)
433  CALL dbcsr_deallocate_matrix_set(matrix_p)
434  CALL dbcsr_deallocate_matrix_set(matrix_h)
435  IF (do_hfx) THEN
436  CALL dbcsr_deallocate_matrix_set(matrix_hfx)
437  END IF
438 
439  CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
440 
441  CALL dbcsr_release_p(fm_deriv)
442  CALL dbcsr_release_p(fm_scaled)
443 
444  DEALLOCATE (occupations)
445  DEALLOCATE (energy_scaling)
446  DEALLOCATE (scaling)
447 
448  CALL dbcsr_release_p(dbcsr_deriv)
449 
450  CALL timestop(handle)
451 
452  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  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  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, rho_ao, tmp_dbcsr
487  TYPE(dbcsr_type), POINTER :: orb_density_matrix, orb_h
488  TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
489  TYPE(pw_c1d_gs_type) :: work_v_gspace
490  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  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  IF (dft_control%sic_method_id .NE. sic_eo) RETURN
503 
504  CALL timeset(routinen, handle)
505 
506  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  virial=virial)
518 
519  CALL qs_rho_get(rho, rho_ao=rho_ao)
520 
521  compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
522  xc_section => section_vals_get_subs_vals(input, "DFT%XC")
523 
524  DO i = 1, SIZE(mo_array) !fm->dbcsr
525  IF (mo_array(i)%use_mo_coeff_b) THEN !fm->dbcsr
526  CALL copy_dbcsr_to_fm(mo_array(i)%mo_coeff_b, &
527  mo_array(i)%mo_coeff) !fm->dbcsr
528  END IF !fm->dbcsr
529  END DO !fm->dbcsr
530 
531  CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
532 
533  ! we have two spins
534  cpassert(SIZE(mo_array, 1) == 2)
535  ! we want uniform occupations
536  CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
537  cpassert(uniform_occupation)
538  CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff, uniform_occupation=uniform_occupation)
539  cpassert(uniform_occupation)
540 
541  NULLIFY (tmp_dbcsr)
542  CALL dbcsr_allocate_matrix_set(tmp_dbcsr, SIZE(mo_derivs, 1))
543  DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
544  !
545  NULLIFY (tmp_dbcsr(i)%matrix)
546  CALL dbcsr_init_p(tmp_dbcsr(i)%matrix)
547  CALL dbcsr_copy(tmp_dbcsr(i)%matrix, mo_derivs(i)%matrix)
548  CALL dbcsr_set(tmp_dbcsr(i)%matrix, 0.0_dp)
549  END DO !fm->dbcsr
550 
551  k_alpha = 0; k_beta = 0
552  SELECT CASE (dft_control%sic_list_id)
553  CASE (sic_list_all)
554 
555  CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
556  CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
557 
558  IF (SIZE(mo_array, 1) > 1) THEN
559  CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
560  CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
561  END IF
562 
563  norb = k_alpha + k_beta
564  ALLOCATE (sic_orbital_list(3, norb))
565 
566  iorb = 0
567  DO i = 1, k_alpha
568  iorb = iorb + 1
569  sic_orbital_list(1, iorb) = 1
570  sic_orbital_list(2, iorb) = i
571  sic_orbital_list(3, iorb) = 1
572  END DO
573  DO i = 1, k_beta
574  iorb = iorb + 1
575  sic_orbital_list(1, iorb) = 2
576  sic_orbital_list(2, iorb) = i
577  IF (SIZE(mo_derivs, 1) == 1) THEN
578  sic_orbital_list(3, iorb) = 1
579  ELSE
580  sic_orbital_list(3, iorb) = 2
581  END IF
582  END DO
583 
584  CASE (sic_list_unpaired)
585  ! we have two spins
586  cpassert(SIZE(mo_array, 1) == 2)
587  ! we have them restricted
588  cpassert(SIZE(mo_derivs, 1) == 1)
589  cpassert(dft_control%restricted)
590 
591  CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
592  CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
593 
594  CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
595  CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
596 
597  norb = k_alpha - k_beta
598  ALLOCATE (sic_orbital_list(3, norb))
599 
600  iorb = 0
601  DO i = k_beta + 1, k_alpha
602  iorb = iorb + 1
603  sic_orbital_list(1, iorb) = 1
604  sic_orbital_list(2, iorb) = i
605  ! we are guaranteed to be restricted
606  sic_orbital_list(3, iorb) = 1
607  END DO
608 
609  CASE DEFAULT
610  cpabort("")
611  END SELECT
612 
613  ! data needed for each of the orbs
614  CALL auxbas_pw_pool%create_pw(orb_rho_r)
615  CALL auxbas_pw_pool%create_pw(tmp_r)
616  CALL auxbas_pw_pool%create_pw(orb_rho_g)
617  CALL auxbas_pw_pool%create_pw(tmp_g)
618  CALL auxbas_pw_pool%create_pw(work_v_gspace)
619  CALL auxbas_pw_pool%create_pw(work_v_rspace)
620 
621  ALLOCATE (orb_density_matrix)
622  CALL dbcsr_copy(orb_density_matrix, rho_ao(1)%matrix, &
623  name="orb_density_matrix")
624  CALL dbcsr_set(orb_density_matrix, 0.0_dp)
625  orb_density_matrix_p%matrix => orb_density_matrix
626 
627  ALLOCATE (orb_h)
628  CALL dbcsr_copy(orb_h, rho_ao(1)%matrix, &
629  name="orb_density_matrix")
630  CALL dbcsr_set(orb_h, 0.0_dp)
631  orb_h_p%matrix => orb_h
632 
633  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  template_fmstruct=mo_coeff%matrix_struct)
637  CALL cp_fm_create(matrix_v, fm_struct_tmp, name="matrix_v")
638  CALL cp_fm_create(matrix_hv, fm_struct_tmp, name="matrix_hv")
639  CALL cp_fm_struct_release(fm_struct_tmp)
640 
641  ALLOCATE (mo_derivs_local(SIZE(mo_array, 1)))
642  DO i = 1, SIZE(mo_array, 1)
643  CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff)
644  CALL cp_fm_create(mo_derivs_local(i), mo_coeff%matrix_struct)
645  END DO
646 
647  ALLOCATE (rho_r(2))
648  rho_r(1) = orb_rho_r
649  rho_r(2) = tmp_r
650  CALL pw_zero(tmp_r)
651 
652  ALLOCATE (rho_g(2))
653  rho_g(1) = orb_rho_g
654  rho_g(2) = tmp_g
655  CALL pw_zero(tmp_g)
656 
657  NULLIFY (vxc)
658  ! now apply to SIC correction to each selected orbital
659  DO iorb = 1, norb
660  ! extract the proper orbital from the mo_coeff
661  CALL get_mo_set(mo_set=mo_array(sic_orbital_list(1, iorb)), mo_coeff=mo_coeff)
662  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  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  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  ks_env=ks_env)
672 
673  ! compute the energy functional for this orbital and its derivative
674 
675  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  energy%hartree = energy%hartree - dft_control%sic_scaling_a*ener
679  IF (.NOT. just_energy) THEN
680  CALL pw_transfer(work_v_gspace, work_v_rspace)
681  CALL pw_scale(work_v_rspace, -dft_control%sic_scaling_a*work_v_rspace%pw_grid%dvol)
682  CALL dbcsr_set(orb_h, 0.0_dp)
683  END IF
684 
685  IF (just_energy) THEN
686  exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
687  pw_pool=xc_pw_pool)
688  ELSE
689  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  pw_pool=xc_pw_pool, compute_virial=compute_virial, virial_xc=virial_xc_tmp)
693  ! add to the existing work_v_rspace
694  CALL pw_axpy(vxc(1), work_v_rspace, -dft_control%sic_scaling_b*vxc(1)%pw_grid%dvol)
695  END IF
696  energy%exc = energy%exc - dft_control%sic_scaling_b*exc
697 
698  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  qs_env=qs_env, calculate_forces=calculate_forces)
702 
703  ! add this to the mo_derivs
704  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  CALL cp_fm_set_all(mo_derivs_local(sic_orbital_list(3, iorb)), 0.0_dp)
707  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  tmp_dbcsr(sic_orbital_list(3, iorb))%matrix)
710  CALL dbcsr_add(mo_derivs(sic_orbital_list(3, iorb))%matrix, &
711  tmp_dbcsr(sic_orbital_list(3, iorb))%matrix, 1.0_dp, 1.0_dp)
712  !
713  ! need to deallocate vxc
714  CALL xc_pw_pool%give_back_pw(vxc(1))
715  CALL xc_pw_pool%give_back_pw(vxc(2))
716  DEALLOCATE (vxc)
717 
718  END IF
719 
720  END DO
721 
722  CALL auxbas_pw_pool%give_back_pw(orb_rho_r)
723  CALL auxbas_pw_pool%give_back_pw(tmp_r)
724  CALL auxbas_pw_pool%give_back_pw(orb_rho_g)
725  CALL auxbas_pw_pool%give_back_pw(tmp_g)
726  CALL auxbas_pw_pool%give_back_pw(work_v_gspace)
727  CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
728 
729  CALL dbcsr_deallocate_matrix(orb_density_matrix)
730  CALL dbcsr_deallocate_matrix(orb_h)
731  CALL cp_fm_release(matrix_v)
732  CALL cp_fm_release(matrix_hv)
733  CALL cp_fm_release(mo_derivs_local)
734  DEALLOCATE (rho_r)
735  DEALLOCATE (rho_g)
736 
737  CALL dbcsr_deallocate_matrix_set(tmp_dbcsr) !fm->dbcsr
738 
739  CALL timestop(handle)
740 
741  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  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  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: store_forces
771  TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
772  TYPE(pw_c1d_gs_type) :: work_rho, work_v
773  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
774  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
775 
776  NULLIFY (mo_array, rho_g)
777 
778  IF (dft_control%sic_method_id == sic_none) RETURN
779  IF (dft_control%sic_method_id == sic_eo) RETURN
780 
781  IF (dft_control%qs_control%gapw) &
782  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  cpassert(dft_control%nspins == 2)
786 
787  CALL auxbas_pw_pool%create_pw(work_rho)
788  CALL auxbas_pw_pool%create_pw(work_v)
789 
790  CALL qs_rho_get(rho, rho_g=rho_g)
791 
792  ! Hartree sic corrections
793  SELECT CASE (dft_control%sic_method_id)
794  CASE (sic_mauri_us, sic_mauri_spz)
795  CALL pw_copy(rho_g(1), work_rho)
796  CALL pw_axpy(rho_g(2), work_rho, alpha=-1._dp)
797  CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
798  CASE (sic_ad)
799  ! find out how many elecs we have
800  CALL get_qs_env(qs_env, mos=mo_array)
801  CALL get_mo_set(mo_set=mo_array(1), nelectron=nelec_a)
802  CALL get_mo_set(mo_set=mo_array(2), nelectron=nelec_b)
803  nelec = nelec_a + nelec_b
804  CALL pw_copy(rho_g(1), work_rho)
805  CALL pw_axpy(rho_g(2), work_rho)
806  scaling = 1.0_dp/real(nelec, kind=dp)
807  CALL pw_scale(work_rho, scaling)
808  CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
809  CASE DEFAULT
810  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  IF (calculate_forces) THEN
816  CALL get_qs_env(qs_env=qs_env, force=force)
817  nforce = 0
818  DO i = 1, SIZE(force)
819  nforce = nforce + SIZE(force(i)%ch_pulay, 2)
820  END DO
821  ALLOCATE (store_forces(3, nforce))
822  nforce = 0
823  DO i = 1, SIZE(force)
824  store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2)) = force(i)%ch_pulay(:, :)
825  force(i)%ch_pulay(:, :) = 0.0_dp
826  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  itype_of_density="SPIN")
836 
837  SELECT CASE (dft_control%sic_method_id)
838  CASE (sic_mauri_us, sic_mauri_spz)
839  full_scaling = -dft_control%sic_scaling_a
840  CASE (sic_ad)
841  full_scaling = -dft_control%sic_scaling_a*nelec
842  CASE DEFAULT
843  cpabort("Unknown sic method id")
844  END SELECT
845  energy%hartree = energy%hartree + full_scaling*ener
846 
847  ! add scaled forces, restoring the old
848  IF (calculate_forces) THEN
849  nforce = 0
850  DO i = 1, SIZE(force)
851  force(i)%ch_pulay(:, :) = force(i)%ch_pulay(:, :)*full_scaling + &
852  store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2))
853  nforce = nforce + SIZE(force(i)%ch_pulay, 2)
854  END DO
855  END IF
856 
857  IF (.NOT. just_energy) THEN
858  ALLOCATE (v_sic_rspace)
859  CALL auxbas_pw_pool%create_pw(v_sic_rspace)
860  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  dft_control%sic_scaling_a*v_sic_rspace%pw_grid%dvol)
864  END IF
865 
866  CALL auxbas_pw_pool%give_back_pw(work_rho)
867  CALL auxbas_pw_pool%give_back_pw(work_v)
868 
869  END SUBROUTINE calc_v_sic_rspace
870 
871 ! **************************************************************************************************
872 !> \brief ...
873 !> \param qs_env ...
874 !> \param rho ...
875 ! **************************************************************************************************
876  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  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r_arr
884  TYPE(cell_type), POINTER :: cell
885  TYPE(cp_logger_type), POINTER :: logger
886  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  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
890  TYPE(section_vals_type), POINTER :: input, scf_section
891 
892  NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, &
893  dft_control, tot_rho_r_arr, rho_ao)
894 
895  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  dft_control=dft_control)
903 
904  CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
905 
906  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  extension=".scfLog")
909 
910  CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao)
911  n_electrons = n_electrons - dft_control%charge
912  tot_rho_r = accurate_sum(tot_rho_r_arr)
913 
914  trace = 0
915  IF (btest(cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES"), cp_p_file)) THEN
916  DO ispin = 1, dft_control%nspins
917  DO img = 1, dft_control%nimages
918  CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp)
919  trace = trace + trace_tmp
920  END DO
921  END DO
922  END IF
923 
924  IF (output_unit > 0) THEN
925  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  "Electronic density on regular grids: ", &
928  tot_rho_r, &
929  tot_rho_r + &
930  REAL(n_electrons, dp), &
931  "Core density on regular grids:", &
932  qs_charges%total_rho_core_rspace, &
933  qs_charges%total_rho_core_rspace - REAL(n_electrons + dft_control%charge, dp)
934  END IF
935  IF (dft_control%qs_control%gapw) THEN
936  tot1_h = qs_charges%total_rho1_hard(1)
937  tot1_s = qs_charges%total_rho1_soft(1)
938  DO ispin = 2, dft_control%nspins
939  tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
940  tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
941  END DO
942  IF (output_unit > 0) THEN
943  WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
944  "Hard and soft densities (Lebedev):", &
945  tot1_h, tot1_s
946  WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
947  "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
948  tot_rho_r + tot1_h - tot1_s, &
949  "Total charge density (r-space): ", &
950  tot_rho_r + tot1_h - tot1_s &
951  + qs_charges%total_rho_core_rspace, &
952  "Total Rho_soft + Rho0_soft (g-space):", &
953  qs_charges%total_rho_gspace
954  END IF
955  qs_charges%background = tot_rho_r + tot1_h - tot1_s + &
956  qs_charges%total_rho_core_rspace
957  ELSE IF (dft_control%qs_control%gapw_xc) THEN
958  tot1_h = qs_charges%total_rho1_hard(1)
959  tot1_s = qs_charges%total_rho1_soft(1)
960  DO ispin = 2, dft_control%nspins
961  tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
962  tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
963  END DO
964  IF (output_unit > 0) THEN
965  WRITE (unit=output_unit, fmt="(/,(T3,A,T41,2F20.10))") &
966  "Hard and soft densities (Lebedev):", &
967  tot1_h, tot1_s
968  WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
969  "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
970  accurate_sum(tot_rho_r_arr) + tot1_h - tot1_s
971  END IF
972  qs_charges%background = tot_rho_r + &
973  qs_charges%total_rho_core_rspace
974  ELSE
975  IF (output_unit > 0) THEN
976  WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
977  "Total charge density on r-space grids: ", &
978  tot_rho_r + &
979  qs_charges%total_rho_core_rspace, &
980  "Total charge density g-space grids: ", &
981  qs_charges%total_rho_gspace
982  END IF
983  qs_charges%background = tot_rho_r + &
984  qs_charges%total_rho_core_rspace
985  END IF
986  IF (output_unit > 0) WRITE (unit=output_unit, fmt="()")
987  qs_charges%background = qs_charges%background/cell%deth
988 
989  CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
990  "PRINT%TOTAL_DENSITIES")
991 
992  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  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  logger => cp_get_default_logger()
1021 
1022  NULLIFY (pw_env)
1023  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1024  psolver = pw_env%poisson_env%parameters%solver
1025 
1026  output_unit = cp_print_key_unit_nr(logger, input, "DFT%SCF%PRINT%DETAILED_ENERGY", &
1027  extension=".scfLog")
1028  IF (output_unit > 0) THEN
1029  IF (dft_control%do_admm) THEN
1030  WRITE (unit=output_unit, fmt="((T3,A,T60,F20.10))") &
1031  "Wfn fit exchange-correlation energy: ", energy%exc_aux_fit
1032  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  "Wfn fit soft/hard atomic rho1 Exc contribution: ", energy%exc1_aux_fit
1035  END IF
1036  END IF
1037  IF (dft_control%do_admm) THEN
1038  IF (psolver .EQ. pw_poisson_implicit) THEN
1039  implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1040  bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1041  SELECT CASE (bc)
1042  CASE (mixed_periodic_bc, mixed_bc)
1043  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1044  "Core Hamiltonian energy: ", energy%core, &
1045  "Hartree energy: ", implicit_ps_ehartree, &
1046  "Electric enthalpy: ", energy%hartree, &
1047  "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  "Core Hamiltonian energy: ", energy%core, &
1051  "Hartree energy: ", energy%hartree, &
1052  "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  "Core Hamiltonian energy: ", energy%core, &
1057  "Hartree energy: ", energy%hartree, &
1058  "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1059  END IF
1060  ELSE
1061 !ZMP to print some variables at each step
1062  IF (dft_control%apply_external_density) THEN
1063  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1064  "DOING ZMP CALCULATION FROM EXTERNAL DENSITY "
1065  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1066  "Core Hamiltonian energy: ", energy%core, &
1067  "Hartree energy: ", energy%hartree
1068  ELSE IF (dft_control%apply_external_vxc) THEN
1069  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1070  "DOING ZMP READING EXTERNAL VXC "
1071  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1072  "Core Hamiltonian energy: ", energy%core, &
1073  "Hartree energy: ", energy%hartree
1074  ELSE
1075  IF (psolver .EQ. pw_poisson_implicit) THEN
1076  implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1077  bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1078  SELECT CASE (bc)
1079  CASE (mixed_periodic_bc, mixed_bc)
1080  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1081  "Core Hamiltonian energy: ", energy%core, &
1082  "Hartree energy: ", implicit_ps_ehartree, &
1083  "Electric enthalpy: ", energy%hartree, &
1084  "Exchange-correlation energy: ", energy%exc
1085  CASE (periodic_bc, neumann_bc)
1086  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1087  "Core Hamiltonian energy: ", energy%core, &
1088  "Hartree energy: ", energy%hartree, &
1089  "Exchange-correlation energy: ", energy%exc
1090  END SELECT
1091  ELSE
1092  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1093  "Core Hamiltonian energy: ", energy%core, &
1094  "Hartree energy: ", energy%hartree, &
1095  "Exchange-correlation energy: ", energy%exc
1096  END IF
1097  END IF
1098  END IF
1099 
1100  IF (dft_control%apply_external_density) THEN
1101  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1102  "Integral of the (density * v_xc): ", energy%exc
1103  END IF
1104 
1105  IF (energy%e_hartree /= 0.0_dp) &
1106  WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1107  "Coulomb (electron-electron) energy: ", energy%e_hartree
1108  IF (energy%dispersion /= 0.0_dp) &
1109  WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1110  "Dispersion energy: ", energy%dispersion
1111  IF (energy%efield /= 0.0_dp) &
1112  WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1113  "Electric field interaction energy: ", energy%efield
1114  IF (energy%gcp /= 0.0_dp) &
1115  WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1116  "gCP energy: ", energy%gcp
1117  IF (dft_control%qs_control%gapw) THEN
1118  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1119  "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit, &
1120  "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
1121  END IF
1122  IF (dft_control%qs_control%gapw_xc) THEN
1123  WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1124  "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit
1125  END IF
1126  IF (dft_control%dft_plus_u) THEN
1127  WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1128  "DFT+U energy:", energy%dft_plus_u
1129  END IF
1130  IF (qs_env%qmmm) THEN
1131  WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1132  "QM/MM Electrostatic energy: ", energy%qmmm_el
1133  IF (qs_env%qmmm_env_qm%image_charge) THEN
1134  WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1135  "QM/MM image charge energy: ", energy%image_charge
1136  END IF
1137  END IF
1138  IF (dft_control%qs_control%mulliken_restraint) THEN
1139  WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1140  "Mulliken restraint (order_p,energy) : ", mulliken_order_p, energy%mulliken
1141  END IF
1142  IF (dft_control%qs_control%ddapc_restraint) THEN
1143  DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
1144  ddapc_order_p = &
1145  dft_control%qs_control%ddapc_restraint_control(n)%ddapc_order_p
1146  WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1147  "DDAPC restraint (order_p,energy) : ", ddapc_order_p, energy%ddapc_restraint(n)
1148  END DO
1149  END IF
1150  IF (dft_control%qs_control%s2_restraint) THEN
1151  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  "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  "DFT%SCF%PRINT%DETAILED_ENERGY")
1159 
1160  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1183  TYPE(dft_control_type), POINTER :: dft_control
1184 
1185  CALL timeset(routinen, handle)
1186 
1187  ! create the matrix using matrix_ks as a template
1188  IF (ASSOCIATED(matrix_vxc)) THEN
1189  CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1190  END IF
1191  CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
1192  ALLOCATE (matrix_vxc(SIZE(matrix_ks)))
1193  DO ispin = 1, SIZE(matrix_ks)
1194  NULLIFY (matrix_vxc(ispin)%matrix)
1195  CALL dbcsr_init_p(matrix_vxc(ispin)%matrix)
1196  CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, &
1197  name="Matrix VXC of spin "//cp_to_string(ispin))
1198  CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp)
1199  END DO
1200 
1201  ! and integrate
1202  CALL get_qs_env(qs_env, dft_control=dft_control)
1203  gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1204  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  gapw=gapw)
1210  ! scale by the volume element... should really become part of integrate_v_rspace
1211  CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw_grid%dvol)
1212  END DO
1213 
1214  CALL timestop(handle)
1215 
1216  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  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  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  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1272  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao, rho_ao_nokp, smat
1273  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, &
1274  matrix_ks_aux_fit_dft, rho_ao_aux, &
1275  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  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  CALL timeset(routinen, handle)
1291 
1292  NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, &
1293  v_rspace, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, &
1294  ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, &
1295  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  vee=vee)
1302 
1303  CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1304  CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
1305  gapw = dft_control%qs_control%gapw
1306  gapw_xc = dft_control%qs_control%gapw_xc
1307  do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
1308 
1309  rigpw = dft_control%qs_control%rigpw
1310  lrigpw = dft_control%qs_control%lrigpw
1311  IF (lrigpw .OR. rigpw) THEN
1312  CALL get_qs_env(qs_env, &
1313  lri_env=lri_env, &
1314  lri_density=lri_density, &
1315  atomic_kind_set=atomic_kind_set)
1316  END IF
1317 
1318  nspins = dft_control%nspins
1319 
1320  ! sum up potentials and integrate
1321  IF (ASSOCIATED(v_rspace_new)) THEN
1322  DO ispin = 1, nspins
1323  IF (gapw_xc) THEN
1324  ! SIC not implemented (or at least not tested)
1325  cpassert(dft_control%sic_method_id == sic_none)
1326  !Only the xc potential, because it has to be integrated with the soft basis
1327  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  rho_ao => rho_ao_kp(ispin, :)
1331  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  gapw=gapw_xc)
1337 
1338  ! Now the Hartree potential to be integrated with the full basis
1339  CALL pw_copy(v_rspace, v_rspace_new(ispin))
1340  ELSE
1341  ! Add v_hartree + v_xc = v_rspace_new
1342  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  IF (dft_control%qs_control%ddapc_explicit_potential) THEN
1345  IF (dft_control%qs_control%ddapc_restraint_is_spin) THEN
1346  IF (ispin == 1) THEN
1347  CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1348  ELSE
1349  CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), -1.0_dp)
1350  END IF
1351  ELSE
1352  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  IF (dft_control%qs_control%cdft) THEN
1357  DO igroup = 1, SIZE(cdft_control%group)
1358  SELECT CASE (cdft_control%group(igroup)%constraint_type)
1359  CASE (cdft_charge_constraint)
1360  csign = 1.0_dp
1361  CASE (cdft_magnetization_constraint)
1362  IF (ispin == 1) THEN
1363  csign = 1.0_dp
1364  ELSE
1365  csign = -1.0_dp
1366  END IF
1367  CASE (cdft_alpha_constraint)
1368  csign = 1.0_dp
1369  IF (ispin == 2) cycle
1370  CASE (cdft_beta_constraint)
1371  csign = 1.0_dp
1372  IF (ispin == 1) cycle
1373  CASE DEFAULT
1374  cpabort("Unknown constraint type.")
1375  END SELECT
1376  CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_new(ispin), &
1377  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  IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
1384  dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1385  CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_new(ispin), dvol)
1386  END IF
1387  ! Add SCCS contribution
1388  IF (dft_control%do_sccs) THEN
1389  CALL pw_axpy(v_sccs_rspace, v_rspace_new(ispin))
1390  END IF
1391  ! External electrostatic potential
1392  IF (dft_control%apply_external_potential) THEN
1393  CALL qmmm_modify_hartree_pot(v_hartree=v_rspace_new(ispin), &
1394  v_qmmm=vee, scale=-1.0_dp)
1395  END IF
1396  IF (do_ppl) THEN
1397  cpassert(.NOT. gapw)
1398  CALL pw_axpy(vppl_rspace, v_rspace_new(ispin), vppl_rspace%pw_grid%dvol)
1399  END IF
1400  ! the electrostatic sic contribution
1401  SELECT CASE (dft_control%sic_method_id)
1402  CASE (sic_none)
1403  !
1404  CASE (sic_mauri_us, sic_mauri_spz)
1405  IF (ispin == 1) THEN
1406  CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1407  ELSE
1408  CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), 1.0_dp)
1409  END IF
1410  CASE (sic_ad)
1411  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  IF (dft_control%apply_embed_pot) THEN
1417  CALL pw_axpy(v_rspace_embed(ispin), v_rspace_new(ispin), v_rspace_embed(ispin)%pw_grid%dvol)
1418  CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1419  END IF
1420  IF (lrigpw) THEN
1421  lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1422  CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1423  DO ikind = 1, nkind
1424  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  lri_v_int, calculate_forces, "LRI_AUX")
1428  DO ikind = 1, nkind
1429  CALL para_env%sum(lri_v_int(ikind)%v_int)
1430  END DO
1431  IF (lri_env%exact_1c_terms) THEN
1432  rho_ao => my_rho(ispin, :)
1433  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  calculate_forces, "ORB")
1437  END IF
1438  IF (lri_env%ppl_ri) THEN
1439  CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1440  END IF
1441  ELSEIF (rigpw) THEN
1442  lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1443  CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1444  DO ikind = 1, nkind
1445  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  lri_v_int, calculate_forces, "RI_HXC")
1449  DO ikind = 1, nkind
1450  CALL para_env%sum(lri_v_int(ikind)%v_int)
1451  END DO
1452  ELSE
1453  rho_ao => my_rho(ispin, :)
1454  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  gapw=gapw)
1460  END IF
1461  CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
1462  END DO ! ispin
1463 
1464  SELECT CASE (dft_control%sic_method_id)
1465  CASE (sic_none)
1466  CASE (sic_mauri_us, sic_mauri_spz, sic_ad)
1467  CALL auxbas_pw_pool%give_back_pw(v_sic_rspace)
1468  DEALLOCATE (v_sic_rspace)
1469  END SELECT
1470  DEALLOCATE (v_rspace_new)
1471 
1472  ELSE
1473  ! not implemented (or at least not tested)
1474  cpassert(dft_control%sic_method_id == sic_none)
1475  cpassert(.NOT. dft_control%qs_control%ddapc_restraint_is_spin)
1476  DO ispin = 1, nspins
1477  ! extra contribution attributed to the dependency of the dielectric constant to the charge density
1478  IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
1479  dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1480  CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace, dvol)
1481  END IF
1482  ! Add SCCS contribution
1483  IF (dft_control%do_sccs) THEN
1484  CALL pw_axpy(v_sccs_rspace, v_rspace)
1485  END IF
1486  ! DFT embedding
1487  IF (dft_control%apply_embed_pot) THEN
1488  CALL pw_axpy(v_rspace_embed(ispin), v_rspace, v_rspace_embed(ispin)%pw_grid%dvol)
1489  CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1490  END IF
1491  IF (lrigpw) THEN
1492  lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1493  CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1494  DO ikind = 1, nkind
1495  lri_v_int(ikind)%v_int = 0.0_dp
1496  END DO
1497  CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
1498  lri_v_int, calculate_forces, "LRI_AUX")
1499  DO ikind = 1, nkind
1500  CALL para_env%sum(lri_v_int(ikind)%v_int)
1501  END DO
1502  IF (lri_env%exact_1c_terms) THEN
1503  rho_ao => my_rho(ispin, :)
1504  ksmat => ks_matrix(ispin, :)
1505  CALL integrate_v_rspace_diagonal(v_rspace, ksmat(1)%matrix, &
1506  rho_ao(1)%matrix, qs_env, &
1507  calculate_forces, "ORB")
1508  END IF
1509  IF (lri_env%ppl_ri) THEN
1510  CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1511  END IF
1512  ELSEIF (rigpw) THEN
1513  lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1514  CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1515  DO ikind = 1, nkind
1516  lri_v_int(ikind)%v_int = 0.0_dp
1517  END DO
1518  CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
1519  lri_v_int, calculate_forces, "RI_HXC")
1520  DO ikind = 1, nkind
1521  CALL para_env%sum(lri_v_int(ikind)%v_int)
1522  END DO
1523  ELSE
1524  rho_ao => my_rho(ispin, :)
1525  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  gapw=gapw)
1532  END IF
1533  END DO
1534  END IF ! ASSOCIATED(v_rspace_new)
1535 
1536  ! **** LRIGPW: KS matrix from integrated potential
1537  IF (lrigpw) THEN
1538  CALL get_qs_env(qs_env, ks_env=ks_env)
1539  CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1540  CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1541  DO ispin = 1, nspins
1542  ksmat => ks_matrix(ispin, :)
1543  CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set, &
1544  cell_to_index=cell_to_index)
1545  END DO
1546  IF (calculate_forces) THEN
1547  CALL calculate_lri_forces(lri_env, lri_density, qs_env, rho_ao_kp, atomic_kind_set)
1548  END IF
1549  ELSEIF (rigpw) THEN
1550  CALL get_qs_env(qs_env, matrix_s=smat)
1551  DO ispin = 1, nspins
1552  CALL calculate_ri_ks_matrix(lri_env, lri_v_int, ks_matrix(ispin, 1)%matrix, &
1553  smat(1)%matrix, atomic_kind_set, ispin)
1554  END DO
1555  IF (calculate_forces) THEN
1556  rho_ao_nokp => rho_ao_kp(:, 1)
1557  CALL calculate_ri_forces(lri_env, lri_density, qs_env, rho_ao_nokp, atomic_kind_set)
1558  END IF
1559  END IF
1560 
1561  IF (ASSOCIATED(v_tau_rspace)) THEN
1562  IF (lrigpw .OR. rigpw) THEN
1563  cpabort("LRIGPW/RIGPW not implemented for meta-GGAs")
1564  END IF
1565  DO ispin = 1, nspins
1566  CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1567 
1568  rho_ao => rho_ao_kp(ispin, :)
1569  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  gapw=gapw)
1575  CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1576  END DO
1577  DEALLOCATE (v_tau_rspace)
1578  END IF
1579 
1580  ! Add contributions from ADMM if requested
1581  IF (dft_control%do_admm) THEN
1582  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  matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft)
1585  CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
1586  IF (ASSOCIATED(v_rspace_new_aux_fit)) THEN
1587  DO ispin = 1, nspins
1588  ! Calculate the xc potential
1589  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  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  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  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  CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
1603  basis_type = "AUX_FIT"
1604  IF (admm_env%do_gapw) THEN
1605  task_list => admm_env%admm_gapw_env%task_list
1606  basis_type = "AUX_FIT_SOFT"
1607  END IF
1608  fadm = 1.0_dp
1609  ! Calculate bare scaling of force according to Merlot, 1. IF: ADMMP, 2. IF: ADMMS,
1610  IF (admm_env%do_admmp) THEN
1611  fadm = admm_env%gsi(ispin)**2
1612  ELSE IF (admm_env%do_admms) THEN
1613  fadm = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)
1614  END IF
1615 
1616  rho_ao => rho_ao_aux(ispin, :)
1617  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  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  DO img = 1, dft_control%nimages
1632  CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, &
1633  matrix_ks_aux_fit(ispin, img)%matrix, 1.0_dp, -1.0_dp)
1634  END DO
1635 
1636  CALL auxbas_pw_pool%give_back_pw(v_rspace_new_aux_fit(ispin))
1637  END DO
1638  DEALLOCATE (v_rspace_new_aux_fit)
1639  END IF
1640  ! Clean up v_tau_rspace_aux_fit, which is actually not needed
1641  IF (ASSOCIATED(v_tau_rspace_aux_fit)) THEN
1642  DO ispin = 1, nspins
1643  CALL auxbas_pw_pool%give_back_pw(v_tau_rspace_aux_fit(ispin))
1644  END DO
1645  DEALLOCATE (v_tau_rspace_aux_fit)
1646  END IF
1647  END IF
1648 
1649  IF (dft_control%apply_embed_pot) DEALLOCATE (v_rspace_embed)
1650 
1651  CALL timestop(handle)
1652 
1653  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  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  REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_ext_r
1682  TYPE(dft_control_type), POINTER :: dft_control
1683  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  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  CALL timeset(routinen, handle)
1695  NULLIFY (auxbas_pw_pool)
1696  NULLIFY (pw_env)
1697  NULLIFY (poisson_env)
1698  NULLIFY (v_rspace_new)
1699  NULLIFY (dft_control)
1700  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  dft_control=dft_control)
1708  CALL pw_env_get(pw_env=pw_env, &
1709  auxbas_pw_pool=auxbas_pw_pool, &
1710  poisson_env=poisson_env)
1711  CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
1712  nspins = 1
1713  ALLOCATE (v_rspace_new(nspins))
1714  CALL auxbas_pw_pool%create_pw(pw=v_rspace_new(1))
1715  CALL auxbas_pw_pool%create_pw(pw=v_xc_rspace)
1716 
1717  CALL pw_zero(v_rspace_new(1))
1718  do_zmp_read = dft_control%apply_external_vxc
1719  IF (do_zmp_read) THEN
1720  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  v_rspace_new(1)%pw_grid%dvol
1723  ELSE
1724  block
1725  REAL(kind=dp) :: factor
1726  TYPE(pw_c1d_gs_type) :: rho_eff_gspace, v_xc_gspace
1727  CALL auxbas_pw_pool%create_pw(pw=rho_eff_gspace)
1728  CALL auxbas_pw_pool%create_pw(pw=v_xc_gspace)
1729  CALL pw_zero(rho_eff_gspace)
1730  CALL pw_zero(v_xc_gspace)
1731  CALL pw_zero(v_xc_rspace)
1732  factor = pw_integrate_function(rho_g(1))
1733  CALL qs_rho_get(qs_env%rho_external, &
1734  rho_g=rho_ext_g, &
1735  tot_rho_r=tot_rho_ext_r)
1736  factor = tot_rho_ext_r(1)/factor
1737 
1738  CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1739  CALL pw_axpy(rho_ext_g(1), rho_eff_gspace, alpha=-1.0_dp)
1740  ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
1741  CALL section_vals_val_get(ext_den_section, "LAMBDA", r_val=lambda)
1742  CALL section_vals_val_get(ext_den_section, "ZMP_CONSTRAINT", i_val=my_val)
1743  CALL section_vals_val_get(ext_den_section, "FERMI_AMALDI", l_val=fermi_amaldi)
1744 
1745  CALL pw_scale(rho_eff_gspace, a=lambda)
1746  nelectron = nelectron_spin(1)
1747  factor = -1.0_dp/nelectron
1748  CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1749 
1750  CALL pw_poisson_solve(poisson_env, rho_eff_gspace, vhartree=v_xc_gspace)
1751  CALL pw_transfer(v_xc_gspace, v_rspace_new(1))
1752  CALL pw_copy(v_rspace_new(1), v_xc_rspace)
1753 
1754  exc = 0.0_dp
1755  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  CALL auxbas_pw_pool%give_back_pw(rho_eff_gspace)
1761  CALL auxbas_pw_pool%give_back_pw(v_xc_gspace)
1762  END block
1763  END IF
1764 
1765  CALL auxbas_pw_pool%give_back_pw(v_xc_rspace)
1766 
1767  CALL timestop(handle)
1768 
1769  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  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  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1796 
1797  CALL timeset(routinen, handle)
1798 
1799  NULLIFY (auxbas_pw_pool)
1800  NULLIFY (pw_env)
1801  NULLIFY (rho_r)
1802  CALL get_qs_env(qs_env=qs_env, &
1803  pw_env=pw_env, &
1804  rho=rho)
1805  CALL pw_env_get(pw_env=pw_env, &
1806  auxbas_pw_pool=auxbas_pw_pool)
1807  CALL qs_rho_get(rho, rho_r=rho_r)
1808  ALLOCATE (v_rspace_embed(dft_control%nspins))
1809 
1810  embed_corr = 0.0_dp
1811 
1812  DO ispin = 1, dft_control%nspins
1813  CALL auxbas_pw_pool%create_pw(pw=v_rspace_embed(ispin))
1814  CALL pw_zero(v_rspace_embed(ispin))
1815 
1816  CALL pw_copy(qs_env%embed_pot, v_rspace_embed(ispin))
1817  embed_corr_local = 0.0_dp
1818 
1819  ! Spin embedding potential in open-shell case
1820  IF (dft_control%nspins .EQ. 2) THEN
1821  IF (ispin .EQ. 1) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), 1.0_dp)
1822  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  embed_corr_local = pw_integral_ab(v_rspace_embed(ispin), rho_r(ispin))
1826 
1827  embed_corr = embed_corr + embed_corr_local
1828 
1829  END DO
1830 
1831  ! If only energy requiested we delete the potential
1832  IF (just_energy) THEN
1833  DO ispin = 1, dft_control%nspins
1834  CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1835  END DO
1836  DEALLOCATE (v_rspace_embed)
1837  END IF
1838 
1839  CALL timestop(handle)
1840 
1841  END SUBROUTINE get_embed_potential_energy
1842 
1843 END MODULE qs_ks_utils
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
admm_types
Types and set/get functions for auxiliary density matrix methods.
Definition: admm_types.F:15
admm_types::get_admm_env
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition: admm_types.F:593
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_sm_fm_multiply
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
Definition: cp_dbcsr_operations.F:744
cp_dbcsr_operations::copy_dbcsr_to_fm
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
Definition: cp_dbcsr_operations.F:208
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
Definition: cp_dbcsr_operations.F:887
cp_dbcsr_operations::copy_fm_to_dbcsr
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Definition: cp_dbcsr_operations.F:118
cp_ddapc
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition: cp_ddapc.F:15
cp_ddapc::cp_ddapc_apply_cd
subroutine, public cp_ddapc_apply_cd(qs_env, rho_tot_gspace, energy, v_hartree_gspace, calculate_forces, Itype_of_density)
Routine to couple/decouple periodic images with the Bloechl scheme.
Definition: cp_ddapc.F:226
cp_fm_struct
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
cp_fm_types::cp_fm_set_all
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition: cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition: cp_output_handling.F:345
hfx_admm_utils
Utilities for hfx and admm methods.
Definition: hfx_admm_utils.F:17
hfx_admm_utils::tddft_hfx_matrix
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data)
Add the hfx contributions to the Hamiltonian.
Definition: hfx_admm_utils.F:2477
hfx_derivatives
Routines to calculate derivatives with respect to basis function origin.
Definition: hfx_derivatives.F:14
hfx_derivatives::derivatives_four_center
subroutine, public derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, irep, use_virial, adiabatic_rescale_factor, resp_only, external_x_data)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Definition: hfx_derivatives.F:117
hfx_types
Types and set/get functions for HFX.
Definition: hfx_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::sic_list_unpaired
integer, parameter, public sic_list_unpaired
Definition: input_constants.F:573
input_constants::sic_mauri_spz
integer, parameter, public sic_mauri_spz
Definition: input_constants.F:568
input_constants::cdft_beta_constraint
integer, parameter, public cdft_beta_constraint
Definition: input_constants.F:766
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition: input_constants.F:766
input_constants::sic_list_all
integer, parameter, public sic_list_all
Definition: input_constants.F:573
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition: input_constants.F:766
input_constants::sic_eo
integer, parameter, public sic_eo
Definition: input_constants.F:568
input_constants::do_ppl_grid
integer, parameter, public do_ppl_grid
Definition: input_constants.F:274
input_constants::do_admm_aux_exch_func_none
integer, parameter, public do_admm_aux_exch_func_none
Definition: input_constants.F:863
input_constants::sic_mauri_us
integer, parameter, public sic_mauri_us
Definition: input_constants.F:568
input_constants::sic_none
integer, parameter, public sic_none
Definition: input_constants.F:568
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition: input_constants.F:766
input_constants::sic_ad
integer, parameter, public sic_ad
Definition: input_constants.F:568
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition: input_section_types.F:1040
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition: kahan_sum.F:29
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition: kinds.F:57
kpoint_types
Types and basic routines needed for a kpoint calculation.
Definition: kpoint_types.F:15
kpoint_types::get_kpoint_info
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Definition: kpoint_types.F:333
lri_environment_methods
Calculates integral matrices for LRIGPW method lri : local resolution of the identity.
Definition: lri_environment_methods.F:18
lri_environment_methods::v_int_ppl_update
subroutine, public v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
...
Definition: lri_environment_methods.F:1189
lri_environment_types
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Definition: lri_environment_types.F:18
lri_forces
Calculates forces for LRIGPW method lri : local resolution of the identity.
Definition: lri_forces.F:16
lri_forces::calculate_lri_forces
subroutine, public calculate_lri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the lri forces
Definition: lri_forces.F:81
lri_forces::calculate_ri_forces
subroutine, public calculate_ri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the ri forces
Definition: lri_forces.F:673
lri_ks_methods
routines that build the Kohn-Sham matrix for the LRIGPW and xc parts
Definition: lri_ks_methods.F:15
lri_ks_methods::calculate_lri_ks_matrix
subroutine, public calculate_lri_ks_matrix(lri_env, lri_v_int, h_matrix, atomic_kind_set, cell_to_index)
update of LRIGPW KS matrix
Definition: lri_ks_methods.F:65
lri_ks_methods::calculate_ri_ks_matrix
subroutine, public calculate_ri_ks_matrix(lri_env, lri_v_int, h_matrix, s_matrix, atomic_kind_set, ispin)
update of RIGPW KS matrix
Definition: lri_ks_methods.F:284
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
ps_implicit_types
Types containing essential information for running implicit (iterative) Poisson solver.
Definition: ps_implicit_types.F:15
ps_implicit_types::neumann_bc
integer, parameter, public neumann_bc
Definition: ps_implicit_types.F:32
ps_implicit_types::mixed_bc
integer, parameter, public mixed_bc
Definition: ps_implicit_types.F:32
ps_implicit_types::mixed_periodic_bc
integer, parameter, public mixed_periodic_bc
Definition: ps_implicit_types.F:32
ps_implicit_types::periodic_bc
integer, parameter, public periodic_bc
Definition: ps_implicit_types.F:32
pw_env_types
container for various plainwaves related things
Definition: pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
pw_methods
Definition: pw_methods.F:25
pw_poisson_methods
Definition: pw_poisson_methods.F:13
pw_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_poisson_types::pw_poisson_implicit
integer, parameter, public pw_poisson_implicit
Definition: pw_poisson_types.F:66
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
pw_types
Definition: pw_types.F:24
qs_cdft_types
Defines CDFT control structures.
Definition: qs_cdft_types.F:14
qs_charges_types
container for information about total charges on the grids
Definition: qs_charges_types.F:14
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition: qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
Definition: qs_collocate_density.F:1710
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition: qs_integrate_potential.F:22
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_ks_qmmm_methods
Definition: qs_ks_qmmm_methods.F:11
qs_ks_qmmm_methods::qmmm_modify_hartree_pot
subroutine, public qmmm_modify_hartree_pot(v_hartree, v_qmmm, scale)
Modify the hartree potential in order to include the QM/MM correction.
Definition: qs_ks_qmmm_methods.F:198
qs_ks_types
Definition: qs_ks_types.F:15
qs_ks_types::get_ks_env
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Definition: qs_ks_types.F:330
qs_ks_utils
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_utils.F:22
qs_ks_utils::print_densities
subroutine, public print_densities(qs_env, rho)
...
Definition: qs_ks_utils.F:877
qs_ks_utils::get_embed_potential_energy
subroutine, public get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, just_energy)
...
Definition: qs_ks_utils.F:1782
qs_ks_utils::low_spin_roks
subroutine, public low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, calculate_forces, auxbas_pw_pool)
do ROKS calculations yielding low spin states
Definition: qs_ks_utils.F:149
qs_ks_utils::sum_up_and_integrate
subroutine, public sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, cdft_control, calculate_forces)
Sum up all potentials defined on the grid and integrate.
Definition: qs_ks_utils.F:1248
qs_ks_utils::print_detailed_energy
subroutine, public print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
Print detailed energies.
Definition: qs_ks_utils.F:1007
qs_ks_utils::calculate_zmp_potential
subroutine, public calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
Calculate the ZMP potential and energy as in Zhao, Morrison Parr PRA 50i, 2138 (1994) V_c^\lambda def...
Definition: qs_ks_utils.F:1669
qs_ks_utils::calc_v_sic_rspace
subroutine, public calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, just_energy, calculate_forces, auxbas_pw_pool)
do sic calculations on the spin density
Definition: qs_ks_utils.F:758
qs_ks_utils::sic_explicit_orbitals
subroutine, public sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, calculate_forces, auxbas_pw_pool)
do sic calculations on explicit orbitals
Definition: qs_ks_utils.F:465
qs_ks_utils::compute_matrix_vxc
subroutine, public compute_matrix_vxc(qs_env, v_rspace, matrix_vxc)
compute matrix_vxc, defined via the potential created by qs_vxc_create ignores things like tau functi...
Definition: qs_ks_utils.F:1174
qs_mo_types
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
qs_mo_types::get_mo_set
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition: qs_rho_types.F:229
task_list_types
types for task lists
Definition: task_list_types.F:14
virial_types
Definition: virial_types.F:13
xc
Exchange and Correlation functional calculations.
Definition: xc.F:17
xc::xc_exc_calc
real(kind=dp) function, public xc_exc_calc(rho_r, rho_g, tau, xc_section, pw_pool)
calculates just the exchange and correlation energy (no vxc)
Definition: xc.F:1449
xc::xc_vxc_pw_create1
subroutine, public xc_vxc_pw_create1(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, xc_section, pw_pool, compute_virial, virial_xc, exc_r)
Exchange and Correlation functional calculations. depending on the selected functional_routine calls ...
Definition: xc.F:150