(git:55939f2)
Loading...
Searching...
No Matches
qs_ks_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief routines that build the Kohn-Sham matrix (i.e calculate the coulomb
10!> and xc parts
11!> \author Fawzi Mohamed
12!> \par History
13!> - 05.2002 moved from qs_scf (see there the history) [fawzi]
14!> - JGH [30.08.02] multi-grid arrays independent from density and potential
15!> - 10.2002 introduced pools, uses updated rho as input,
16!> removed most temporary variables, renamed may vars,
17!> began conversion to LSD [fawzi]
18!> - 10.2004 moved calculate_w_matrix here [Joost VandeVondele]
19!> introduced energy derivative wrt MOs [Joost VandeVondele]
20!> - SCCS implementation (16.10.2013,MK)
21! **************************************************************************************************
22MODULE qs_ks_methods
23 USE accint_weights_forces, ONLY: accint_weight_force
24 USE admm_dm_methods, ONLY: admm_dm_calc_rho_aux,&
25 admm_dm_merge_ks_matrix
26 USE admm_methods, ONLY: admm_mo_calc_rho_aux,&
27 admm_mo_calc_rho_aux_kp,&
28 admm_mo_merge_ks_matrix,&
29 admm_update_ks_atom,&
30 calc_admm_mo_derivatives,&
31 calc_admm_ovlp_forces,&
32 calc_admm_ovlp_forces_kp
33 USE admm_types, ONLY: admm_type,&
34 get_admm_env
35 USE atomic_kind_types, ONLY: atomic_kind_type,&
36 get_atomic_kind_set
37 USE cell_types, ONLY: cell_type
38 USE cp_control_types, ONLY: dft_control_type
39 USE cp_dbcsr_api, ONLY: &
40 dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_filter, dbcsr_get_info, dbcsr_multiply, &
41 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
42 dbcsr_type_symmetric
43 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
44 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
45 dbcsr_copy_columns_hack
46 USE cp_ddapc, ONLY: qs_ks_ddapc
47 USE cp_fm_types, ONLY: cp_fm_type
48 USE cp_log_handling, ONLY: cp_get_default_logger,&
49 cp_logger_get_default_io_unit,&
50 cp_logger_type
51 USE cp_output_handling, ONLY: cp_p_file,&
52 cp_print_key_should_output
53 USE dft_plus_u, ONLY: plus_u
54 USE gapw_gspace_reconstruction, ONLY: integrate_gapw_composite_vxc_gspace
55 USE gce_methods, ONLY: planar_averaged_v_hartree_3d,&
56 planar_counter_charge
57 USE hartree_local_methods, ONLY: vh_1c_gg_integrals
58 USE hartree_local_types, ONLY: ecoul_1center_type
59 USE hfx_ace_methods, ONLY: hfx_ace_ks_matrix
60 USE hfx_admm_utils, ONLY: hfx_admm_init,&
61 hfx_ks_matrix,&
62 hfx_ks_matrix_kp
63 USE input_constants, ONLY: do_ppl_grid,&
64 outer_scf_becke_constraint,&
65 outer_scf_hirshfeld_constraint,&
66 smeagol_runtype_emtransport
67 USE input_section_types, ONLY: section_vals_get,&
68 section_vals_get_subs_vals,&
69 section_vals_type,&
70 section_vals_val_get
71 USE kg_correction, ONLY: kg_ekin_subset
72 USE kinds, ONLY: default_string_length,&
73 dp
74 USE kpoint_types, ONLY: get_kpoint_info,&
75 kpoint_type
76 USE lri_environment_methods, ONLY: v_int_ppl_energy
77 USE lri_environment_types, ONLY: lri_density_type,&
78 lri_environment_type,&
79 lri_kind_type
80 USE mathlib, ONLY: abnormal_value
81 USE message_passing, ONLY: mp_para_env_type
82 USE particle_types, ONLY: particle_type
83 USE pw_env_types, ONLY: pw_env_get,&
84 pw_env_type
85 USE pw_methods, ONLY: pw_axpy,&
86 pw_copy,&
87 pw_integral_ab,&
88 pw_integrate_function,&
89 pw_scale,&
90 pw_transfer,&
91 pw_zero
92 USE pw_poisson_methods, ONLY: pw_poisson_solve
93 USE pw_poisson_types, ONLY: pw_poisson_implicit,&
94 pw_poisson_type
95 USE pw_pool_types, ONLY: pw_pool_type
96 USE pw_types, ONLY: pw_c1d_gs_type,&
97 pw_r3d_rs_type
98 USE qmmm_image_charge, ONLY: add_image_pot_to_hartree_pot,&
99 calculate_image_pot,&
100 integrate_potential_devga_rspace
101 USE qs_cdft_types, ONLY: cdft_control_type
102 USE qs_cdft_utils, ONLY: cdft_constraint_print
103 USE qs_charges_types, ONLY: qs_charges_type
104 USE qs_core_energies, ONLY: calculate_ptrace
105 USE qs_dftb_matrices, ONLY: build_dftb_ks_matrix
106 USE qs_efield_berry, ONLY: qs_efield_berry_phase
107 USE qs_efield_local, ONLY: qs_efield_local_operator
108 USE qs_energy_types, ONLY: qs_energy_type
109 USE qs_environment_types, ONLY: get_qs_env,&
110 qs_environment_type
111 USE qs_force_types, ONLY: qs_force_type
112 USE qs_gapw_densities, ONLY: prepare_gapw_den
113 USE qs_harris_types, ONLY: harris_type
114 USE qs_harris_utils, ONLY: harris_set_potentials
115 USE qs_integrate_potential, ONLY: integrate_ppl_rspace,&
116 integrate_rho_nlcc,&
117 integrate_v_core_rspace,&
118 integrate_v_rspace
119 USE qs_kind_types, ONLY: qs_kind_type
120 USE qs_ks_apply_restraints, ONLY: qs_ks_cdft_constraint,&
121 qs_ks_mulliken_restraint,&
122 qs_ks_s2_restraint
123 USE qs_ks_atom, ONLY: update_ks_atom
124 USE qs_ks_qmmm_methods, ONLY: qmmm_calculate_energy,&
125 qmmm_modify_hartree_pot
126 USE qs_ks_types, ONLY: qs_ks_env_type,&
127 set_ks_env
128 USE qs_ks_utils, ONLY: &
129 calc_v_sic_rspace, calculate_zmp_potential, check_sum_energies, compute_matrix_vxc, &
130 compute_matrix_vxc_kp, get_embed_potential_energy, low_spin_roks, print_densities, &
131 print_detailed_energy, sic_explicit_orbitals, sum_up_and_integrate
132 USE qs_local_rho_types, ONLY: local_rho_type
133 USE qs_mo_types, ONLY: get_mo_set,&
134 mo_set_type
135 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
136 USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
137 USE qs_rho_types, ONLY: qs_rho_get,&
138 qs_rho_type
139 USE qs_sccs, ONLY: sccs,&
140 sccs_core_fill_coefficients
141 USE qs_scf_types, ONLY: qs_scf_env_type
142 USE qs_vxc, ONLY: qs_vxc_create
143 USE qs_vxc_atom, ONLY: calculate_vxc_atom,&
144 gapw_cdft_one_center
145 USE rtp_admm_methods, ONLY: rtp_admm_calc_rho_aux,&
146 rtp_admm_merge_ks_matrix
147 USE se_fock_matrix, ONLY: build_se_fock_matrix
148 USE skala_gpw_functional, ONLY: &
149 ensure_native_skala_grid_scope, get_gauxc_section, &
150 native_skala_gapw_atom_composite_requested, native_skala_gapw_composite_direct_ao, &
151 native_skala_gapw_composite_reference, native_skala_uses_atom_composite_grid, &
152 skala_gapw_representation, xc_section_uses_gauxc_model, xc_section_uses_native_skala_grid
153 USE smeagol_interface, ONLY: smeagol_shift_v_hartree
154 USE string_utilities, ONLY: uppercase
155 USE surface_dipole, ONLY: calc_dipsurf_potential
156 USE tblite_ks_matrix, ONLY: build_tblite_ks_matrix
157 USE virial_types, ONLY: virial_type
158 USE xc_gauxc_functional, ONLY: apply_gauxc,&
159 gauxc_gapw_all_all_electron,&
160 gauxc_gapw_all_pseudopotentials,&
161 gauxc_gapw_has_paw_pseudopotentials
162 USE xc_input_constants, ONLY: skala_gapw_direct_valence,&
163 skala_gapw_paw_one_center
164 USE xtb_ks_matrix, ONLY: build_xtb_ks_matrix
165#include "./base/base_uses.f90"
166
167 IMPLICIT NONE
168
169 PRIVATE
170
171 LOGICAL, PARAMETER :: debug_this_module = .true.
172 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_methods'
173
174 PUBLIC :: calc_rho_tot_gspace, qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, &
175 qs_ks_allocate_basics, evaluate_core_matrix_traces, rebuild_ks_matrix
176
177CONTAINS
178
179! **************************************************************************************************
180!> \brief routine where the real calculations are made: the
181!> KS matrix is calculated
182!> \param qs_env the qs_env to update
183!> \param calculate_forces if true calculate the quantities needed
184!> to calculate the forces. Defaults to false.
185!> \param just_energy if true updates the energies but not the
186!> ks matrix. Defaults to false
187!> \param print_active ...
188!> \param ext_ks_matrix ...
189!> \param ext_xc_section ...
190!> \par History
191!> 06.2002 moved from qs_scf to qs_ks_methods, use of ks_env
192!> new did_change scheme [fawzi]
193!> 10.2002 introduced pools, uses updated rho as input, LSD [fawzi]
194!> 10.2004 build_kohn_sham matrix now also computes the derivatives
195!> of the total energy wrt to the MO coefs, if instructed to
196!> do so. This appears useful for orbital dependent functionals
197!> where the KS matrix alone (however this might be defined)
198!> does not contain the info to construct this derivative.
199!> \author Matthias Krack
200!> \note
201!> make rho, energy and qs_charges optional, defaulting
202!> to qs_env components?
203! **************************************************************************************************
204 SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
205 print_active, ext_ks_matrix, ext_xc_section)
206 TYPE(qs_environment_type), POINTER :: qs_env
207 LOGICAL, INTENT(in) :: calculate_forces, just_energy
208 LOGICAL, INTENT(IN), OPTIONAL :: print_active
209 TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
210 POINTER :: ext_ks_matrix
211 TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
212
213 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_ks_build_kohn_sham_matrix'
214
215 CHARACTER(len=default_string_length) :: gauxc_model_name, name
216 INTEGER :: ace_rebuild_frequency, atom_a, gapw_representation, handle, iatom, ikind, img, &
217 ispin, natom, nimages, nspins, output_unit
218 INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
219 LOGICAL :: ace_active, do_adiabatic_rescaling, do_ddapc, do_hfx, do_kpoints, do_ppl, dokp, &
220 gapw, gapw_all_all_electron, gapw_all_pseudopotentials, gapw_mixed_core_representations, &
221 gapw_xc, gauxc_model_none, just_energy_xc, lrigpw, my_print, &
222 native_direct_valence_atom_composite, native_gapw_all_electron_atom_composite, &
223 native_gapw_atom_composite_reference, native_gapw_atom_composite_requested, &
224 native_gapw_composite_direct_ao, native_gapw_composite_reference, &
225 native_gapw_periodic_atom_composite, native_grid_diagnostics, native_grid_use_cuda, &
226 native_skala_evaluator, native_skala_restore_exc
227 LOGICAL :: paw_composite_representation, rigpw, use_gauxc_matrix, use_virial
228 LOGICAL, SAVE :: native_grid_cpu_kpoints_warned = .false.
229 REAL(kind=dp) :: ecore_ppl, edisp, ee_ener, ekin_mol, mulliken_order_p, &
230 native_skala_composite_contraction, native_skala_composite_rho_contraction, &
231 native_skala_composite_tau_contraction, native_skala_exc_scf, native_skala_total_scf, &
232 vscale
233 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cdft_gapw_values
234 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cdft_charge_correction, &
235 cdft_electronic_charge, &
236 native_skala_atom_force
237 REAL(kind=dp), DIMENSION(3, 3) :: h_stress, native_skala_composite_virial, &
238 pv_loc
239 TYPE(admm_type), POINTER :: admm_env
240 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
241 TYPE(cdft_control_type), POINTER :: cdft_control
242 TYPE(cell_type), POINTER :: cell
243 TYPE(cp_logger_type), POINTER :: logger
244 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_vhxc, matrix_vxc, mo_derivs
245 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, ks_matrix_im, matrix_h, &
246 matrix_h_im, matrix_s, matrix_vxc_kp, &
247 my_rho, rho_ao
248 TYPE(dft_control_type), POINTER :: dft_control
249 TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
250 TYPE(harris_type), POINTER :: harris_env
251 TYPE(kpoint_type), POINTER :: kpoints
252 TYPE(local_rho_type), POINTER :: local_rho_set
253 TYPE(lri_density_type), POINTER :: lri_density
254 TYPE(lri_environment_type), POINTER :: lri_env
255 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
256 TYPE(mp_para_env_type), POINTER :: para_env
257 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
258 POINTER :: sab_orb
259 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
260 TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
261 TYPE(pw_c1d_gs_type), POINTER :: rho_core
262 TYPE(pw_env_type), POINTER :: pw_env
263 TYPE(pw_poisson_type), POINTER :: poisson_env
264 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
265 TYPE(pw_r3d_rs_type), DIMENSION(2) :: v_hxc_rspace
266 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace_embed, v_rspace_new, &
267 v_rspace_new_aux_fit, v_tau_rspace, &
268 v_tau_rspace_aux_fit
269 TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs, rho_nlcc, rhoz_cneo_s_rs, v_hartree_rspace, &
270 v_sccs_rspace, v_sic_rspace, v_spin_ddapc_rest_r, vee, vppl_rspace
271 TYPE(qs_energy_type), POINTER :: energy
272 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
273 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
274 TYPE(qs_ks_env_type), POINTER :: ks_env
275 TYPE(qs_rho_type), POINTER :: rho, rho1, rho_struct, rho_xc
276 TYPE(section_vals_type), POINTER :: ace_section, &
277 adiabatic_rescaling_section, &
278 gauxc_section, hfx_sections, input, &
279 scf_section, xc_section
280 TYPE(virial_type), POINTER :: virial
281
282 CALL timeset(routinen, handle)
283 NULLIFY (admm_env, atomic_kind_set, cell, dft_control, force, logger, mo_derivs, my_rho, &
284 rho_struct, para_env, pw_env, virial, vppl_rspace, &
285 ace_section, &
286 adiabatic_rescaling_section, hfx_sections, input, scf_section, &
287 xc_section, gauxc_section, matrix_h, matrix_h_im, matrix_s, auxbas_pw_pool, poisson_env, &
288 v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, matrix_vxc, &
289 matrix_vxc_kp, matrix_vhxc, &
290 vee, rho_nlcc, ks_env, ks_matrix, ks_matrix_im, rho, energy, rho_xc, rho_r, rho_ao, &
291 rho_core, sab_orb, particle_set, qs_kind_set, kpoints)
292
293 cpassert(ASSOCIATED(qs_env))
294
295 logger => cp_get_default_logger()
296 my_print = .true.
297 IF (PRESENT(print_active)) my_print = print_active
298 use_gauxc_matrix = .false.
299 native_gapw_atom_composite_reference = .false.
300 native_gapw_atom_composite_requested = .false.
301 native_gapw_all_electron_atom_composite = .false.
302 native_direct_valence_atom_composite = .false.
303 gapw_all_all_electron = .false.
304 gapw_all_pseudopotentials = .false.
305 gapw_mixed_core_representations = .false.
306 native_gapw_composite_direct_ao = .false.
307 native_gapw_composite_reference = .false.
308 native_gapw_periodic_atom_composite = .false.
309 native_grid_diagnostics = .false.
310 native_skala_composite_contraction = 0.0_dp
311 native_skala_composite_rho_contraction = 0.0_dp
312 native_skala_composite_tau_contraction = 0.0_dp
313 native_skala_restore_exc = .false.
314
315 CALL get_qs_env(qs_env, &
316 ks_env=ks_env, &
317 dft_control=dft_control, &
318 matrix_h_kp=matrix_h, &
319 matrix_h_im_kp=matrix_h_im, &
320 matrix_s_kp=matrix_s, &
321 matrix_ks_kp=ks_matrix, &
322 matrix_ks_im_kp=ks_matrix_im, &
323 matrix_vxc=matrix_vxc, &
324 matrix_vhxc=matrix_vhxc, &
325 matrix_vxc_kp=matrix_vxc_kp, &
326 pw_env=pw_env, &
327 cell=cell, &
328 atomic_kind_set=atomic_kind_set, &
329 para_env=para_env, &
330 input=input, &
331 virial=virial, &
332 v_hartree_rspace=v_hartree_rspace, &
333 vee=vee, &
334 rho_nlcc=rho_nlcc, &
335 rho=rho, &
336 rho_core=rho_core, &
337 rho_xc=rho_xc, &
338 sab_orb=sab_orb, &
339 energy=energy, &
340 force=force, &
341 kpoints=kpoints, &
342 do_kpoints=do_kpoints, &
343 particle_set=particle_set, &
344 qs_kind_set=qs_kind_set, &
345 natom=natom)
346
347 ! Reject force/stress requests before assembling even the baseline KS matrix.
348 IF (dft_control%low_spin_roks .AND. dft_control%qs_control%gapw) THEN
349 IF (.NOT. virial%pv_numer .AND. &
350 (virial%pv_calculate .OR. (virial%pv_availability .AND. calculate_forces))) THEN
351 CALL cp_abort(__location__, "GAPW low spin ROKS stress is not implemented.")
352 END IF
353 IF (calculate_forces) THEN
354 CALL cp_abort(__location__, "GAPW low spin ROKS forces are not implemented.")
355 END IF
356 END IF
357
358 CALL qs_rho_get(rho, rho_r=rho_r, rho_ao_kp=rho_ao)
359
360 nimages = dft_control%nimages
361 nspins = dft_control%nspins
362
363 ! remap pointer to allow for non-kpoint external ks matrix
364 IF (PRESENT(ext_ks_matrix)) ks_matrix(1:nspins, 1:1) => ext_ks_matrix(1:nspins)
365
366 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
367
368 adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
369 CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
370 just_energy_xc = just_energy
371 IF (do_adiabatic_rescaling) THEN
372 !! If we perform adiabatic rescaling, the xc potential has to be scaled by the xc- and
373 !! HFX-energy. Thus, let us first calculate the energy
374 just_energy_xc = .true.
375 END IF
376
377 cpassert(ASSOCIATED(matrix_h))
378 cpassert(ASSOCIATED(matrix_s))
379 cpassert(ASSOCIATED(rho))
380 cpassert(ASSOCIATED(pw_env))
381 cpassert(SIZE(ks_matrix, 1) > 0)
382 dokp = (nimages > 1)
383
384 ! Setup the possible usage of DDAPC charges
385 do_ddapc = dft_control%qs_control%ddapc_restraint .OR. &
386 qs_env%cp_ddapc_ewald%do_decoupling .OR. &
387 qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
388 qs_env%cp_ddapc_ewald%do_solvation
389
390 ! Check if LRIGPW is used
391 lrigpw = dft_control%qs_control%lrigpw
392 rigpw = dft_control%qs_control%rigpw
393 IF (rigpw) THEN
394 cpassert(nimages == 1)
395 END IF
396 IF (lrigpw .AND. rigpw) THEN
397 cpabort(" LRI and RI are not compatible")
398 END IF
399
400 ! Check for GAPW method : additional terms for local densities
401 gapw = dft_control%qs_control%gapw
402 gapw_xc = dft_control%qs_control%gapw_xc
403 IF (gapw_xc .AND. gapw) THEN
404 cpabort(" GAPW and GAPW_XC are not compatible")
405 END IF
406 IF ((gapw .AND. lrigpw) .OR. (gapw_xc .AND. lrigpw)) THEN
407 cpabort(" GAPW/GAPW_XC and LRIGPW are not compatible")
408 END IF
409 IF ((gapw .AND. rigpw) .OR. (gapw_xc .AND. rigpw)) THEN
410 cpabort(" GAPW/GAPW_XC and RIGPW are not compatible")
411 END IF
412
413 do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
414 IF (do_ppl) THEN
415 cpassert(.NOT. gapw)
416 CALL get_qs_env(qs_env=qs_env, vppl=vppl_rspace)
417 END IF
418
419 IF (gapw_xc) THEN
420 cpassert(ASSOCIATED(rho_xc))
421 END IF
422
423 ! gets the tmp grids
424 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
425
426 IF (gapw .AND. (poisson_env%parameters%solver == pw_poisson_implicit)) THEN
427 cpabort("The implicit Poisson solver cannot be used in conjunction with GAPW.")
428 END IF
429
430 ! *** Prepare densities for gapw ***
431 IF (gapw .OR. gapw_xc) THEN
432 CALL prepare_gapw_den(qs_env, do_rho0=(.NOT. gapw_xc))
433 END IF
434
435 ! Calculate the Hartree potential
436 CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
437 CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
438
439 scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
440 IF (btest(cp_print_key_should_output(logger%iter_info, scf_section, &
441 "PRINT%DETAILED_ENERGY"), &
442 cp_p_file) .AND. &
443 (.NOT. gapw) .AND. (.NOT. gapw_xc) .AND. &
444 (.NOT. (poisson_env%parameters%solver == pw_poisson_implicit))) THEN
445 CALL pw_zero(rho_tot_gspace)
446 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density=.true.)
447 CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%e_hartree, &
448 v_hartree_gspace)
449 CALL pw_zero(rho_tot_gspace)
450 CALL pw_zero(v_hartree_gspace)
451 END IF
452
453 ! Get the total density in g-space [ions + electrons]
454 CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
455
456 IF (qs_env%scf_control%gce%do_gce .AND. .NOT. dft_control%do_pcc) THEN
457 cpabort("GCE requires DFT%PLANAR_COUNTER_CHARGE to define the countercharge plane.")
458 END IF
459
460 ! Add the planar counter charge density
461 IF (dft_control%do_pcc) THEN
462 CALL planar_counter_charge(rho_tot_gspace, dft_control%pcc_control, auxbas_pw_pool)
463 END IF
464
465 IF (my_print) THEN
466 CALL print_densities(qs_env, rho)
467 END IF
468
469 IF (dft_control%do_sccs) THEN
470 ! Self-consistent continuum solvation (SCCS) model
471 NULLIFY (v_sccs_rspace)
472 ALLOCATE (v_sccs_rspace)
473 CALL auxbas_pw_pool%create_pw(v_sccs_rspace)
474
475 IF (poisson_env%parameters%solver == pw_poisson_implicit) THEN
476 cpabort("The implicit Poisson solver cannot be used together with SCCS.")
477 END IF
478
479 IF (use_virial .AND. calculate_forces) THEN
480 CALL sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs_rspace, &
481 h_stress=h_stress)
482 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
483 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
484 ELSE
485 CALL sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs_rspace)
486 END IF
487 ELSE
488 ! Getting the Hartree energy and Hartree potential. Also getting the stress tensor
489 ! from the Hartree term if needed. No nuclear force information here
490 IF (use_virial .AND. calculate_forces) THEN
491 h_stress(:, :) = 0.0_dp
492 CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%hartree, &
493 v_hartree_gspace, h_stress=h_stress, &
494 rho_core=rho_core)
495 virial%pv_ehartree = virial%pv_ehartree + h_stress/real(para_env%num_pe, dp)
496 virial%pv_virial = virial%pv_virial + h_stress/real(para_env%num_pe, dp)
497 ELSE
498 CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%hartree, &
499 v_hartree_gspace, rho_core=rho_core)
500 END IF
501 END IF
502
503 IF (dft_control%do_paep .OR. qs_env%scf_control%gce%do_gce) THEN
504 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
505 CALL planar_averaged_v_hartree_3d(v_hartree_rspace, dft_control, qs_env%scf_control%gce%do_gce, &
506 qs_env%scf_control%gce%ref_esp, para_env)
507 END IF
508
509 ! In case decouple periodic images and/or apply restraints to charges
510 IF (do_ddapc) THEN
511 CALL qs_ks_ddapc(qs_env, auxbas_pw_pool, rho_tot_gspace, v_hartree_gspace, &
512 v_spin_ddapc_rest_r, energy, calculate_forces, ks_matrix, &
513 just_energy)
514 ELSE
515 dft_control%qs_control%ddapc_explicit_potential = .false.
516 dft_control%qs_control%ddapc_restraint_is_spin = .false.
517 IF (.NOT. just_energy) THEN
518 CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
519 CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
520 END IF
521 END IF
522 CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
523
524 IF (dft_control%correct_surf_dip) THEN
525 IF (dft_control%surf_dip_correct_switch) THEN
526 CALL calc_dipsurf_potential(qs_env, energy)
527 energy%hartree = energy%hartree + energy%surf_dipole
528 END IF
529 END IF
530
531 ! SIC
532 CALL calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, &
533 just_energy, calculate_forces, auxbas_pw_pool)
534
535 ! Check if CDFT constraint is needed
536 CALL qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
537
538 ! Adds the External Potential if requested
539 IF (dft_control%apply_external_potential) THEN
540 ! Compute the energy due to the external potential
541 ee_ener = 0.0_dp
542 DO ispin = 1, nspins
543 ee_ener = ee_ener + pw_integral_ab(rho_r(ispin), vee)
544 END DO
545 IF (.NOT. just_energy) THEN
546 IF (gapw) THEN
547 CALL get_qs_env(qs_env=qs_env, &
548 rho0_s_rs=rho0_s_rs, &
549 rhoz_cneo_s_rs=rhoz_cneo_s_rs)
550 cpassert(ASSOCIATED(rho0_s_rs))
551 IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
552 CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs)
553 END IF
554 ee_ener = ee_ener + pw_integral_ab(rho0_s_rs, vee)
555 IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
556 CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs, -1.0_dp)
557 END IF
558 END IF
559 END IF
560 ! the sign accounts for the charge of the electrons
561 energy%ee = -ee_ener
562 END IF
563
564 ! Adds the QM/MM potential
565 IF (qs_env%qmmm) THEN
566 CALL qmmm_calculate_energy(qs_env=qs_env, &
567 rho=rho_r, &
568 v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, &
569 qmmm_energy=energy%qmmm_el)
570 IF (qs_env%qmmm_env_qm%image_charge) THEN
571 CALL calculate_image_pot(v_hartree_rspace=v_hartree_rspace, &
572 rho_hartree_gspace=rho_tot_gspace, &
573 energy=energy, &
574 qmmm_env=qs_env%qmmm_env_qm, &
575 qs_env=qs_env)
576 IF (.NOT. just_energy) THEN
577 CALL add_image_pot_to_hartree_pot(v_hartree=v_hartree_rspace, &
578 v_metal=qs_env%ks_qmmm_env%v_metal_rspace, &
579 qs_env=qs_env)
580 IF (calculate_forces) THEN
581 CALL integrate_potential_devga_rspace( &
582 potential=v_hartree_rspace, coeff=qs_env%image_coeff, &
583 forces=qs_env%qmmm_env_qm%image_charge_pot%image_forcesMM, &
584 qmmm_env=qs_env%qmmm_env_qm, qs_env=qs_env)
585 END IF
586 END IF
587 CALL qs_env%ks_qmmm_env%v_metal_rspace%release()
588 DEALLOCATE (qs_env%ks_qmmm_env%v_metal_rspace)
589 END IF
590 IF (.NOT. just_energy) THEN
591 CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, &
592 v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, scale=1.0_dp)
593 END IF
594 END IF
595 CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
596
597 ! SMEAGOL interface
598 IF (dft_control%smeagol_control%smeagol_enabled .AND. &
599 dft_control%smeagol_control%run_type == smeagol_runtype_emtransport) THEN
600 cpassert(ASSOCIATED(dft_control%smeagol_control%aux))
601 CALL smeagol_shift_v_hartree(v_hartree_rspace, cell, &
602 dft_control%smeagol_control%aux%HartreeLeadsLeft, &
603 dft_control%smeagol_control%aux%HartreeLeadsRight, &
604 dft_control%smeagol_control%aux%HartreeLeadsBottom, &
605 dft_control%smeagol_control%aux%VBias, &
606 dft_control%smeagol_control%aux%minL, &
607 dft_control%smeagol_control%aux%maxR, &
608 dft_control%smeagol_control%aux%isexplicit_maxR, &
609 dft_control%smeagol_control%aux%isexplicit_HartreeLeadsBottom)
610 END IF
611
612 ! calculate the density matrix for the fitted mo_coeffs
613 IF (dft_control%do_admm) THEN
614 IF (PRESENT(ext_xc_section)) THEN
615 CALL hfx_admm_init(qs_env, calculate_forces, ext_xc_section)
616 ELSE
617 CALL hfx_admm_init(qs_env, calculate_forces)
618 END IF
619
620 IF (dft_control%do_admm_mo) THEN
621 IF (qs_env%run_rtp) THEN
622 CALL rtp_admm_calc_rho_aux(qs_env)
623 ELSE
624 IF (dokp) THEN
625 CALL admm_mo_calc_rho_aux_kp(qs_env)
626 ELSE
627 CALL admm_mo_calc_rho_aux(qs_env)
628 END IF
629 END IF
630 ELSE IF (dft_control%do_admm_dm) THEN
631 CALL admm_dm_calc_rho_aux(qs_env)
632 END IF
633 END IF
634
635 ! only activate stress calculation if
636 IF (use_virial .AND. calculate_forces) virial%pv_calculate = .true.
637
638 ! *** calculate the xc potential on the pw density ***
639 ! *** associates v_rspace_new if the xc potential needs to be computed.
640 ! If we do wavefunction fitting, we need the vxc_potential in the auxiliary basis set
641 IF (dft_control%do_admm) THEN
642 CALL get_qs_env(qs_env, admm_env=admm_env)
643 xc_section => admm_env%xc_section_aux
644 CALL get_admm_env(admm_env, rho_aux_fit=rho_struct)
645
646 ! here we ignore a possible vdW section in admm_env%xc_section_aux
647 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
648 vxc_rho=v_rspace_new_aux_fit, vxc_tau=v_tau_rspace_aux_fit, exc=energy%exc_aux_fit, &
649 just_energy=just_energy_xc)
650
651 IF (admm_env%do_gapw) THEN
652 !compute the potential due to atomic densities
653 CALL calculate_vxc_atom(qs_env, energy_only=just_energy_xc, exc1=energy%exc1_aux_fit, &
654 kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
655 xc_section_external=xc_section, &
656 rho_atom_set_external=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
657 calculate_forces=calculate_forces)
658
659 END IF
660
661 NULLIFY (rho_struct)
662
663 IF (use_virial .AND. calculate_forces) THEN
664 vscale = 1.0_dp
665 !Note: ADMMS and ADMMP stress tensor only for closed-shell calculations
666 IF (admm_env%do_admms) vscale = admm_env%gsi(1)**(2.0_dp/3.0_dp)
667 IF (admm_env%do_admmp) vscale = admm_env%gsi(1)**2
668 virial%pv_exc = virial%pv_exc - vscale*virial%pv_xc
669 virial%pv_virial = virial%pv_virial - vscale*virial%pv_xc
670 ! virial%pv_xc will be zeroed in the xc routines
671 END IF
672 xc_section => admm_env%xc_section_primary
673 ELSE
674 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
675 ! build ks matrix with an xc section potentially different from the one defined in input
676 IF (PRESENT(ext_xc_section)) xc_section => ext_xc_section
677 END IF
678
679 IF (gapw_xc) THEN
680 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
681 ELSE
682 CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
683 END IF
684
685 ! zmp
686 IF (dft_control%apply_external_density .OR. dft_control%apply_external_vxc) THEN
687 energy%exc = 0.0_dp
688 CALL calculate_zmp_potential(qs_env, v_rspace_new, rho, exc=energy%exc)
689 ELSE
690 ! Embedding potential (runs regardless of XC method)
691 IF (dft_control%apply_embed_pot) THEN
692 NULLIFY (v_rspace_embed)
693 energy%embed_corr = 0.0_dp
694 CALL get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, &
695 energy%embed_corr, just_energy)
696 END IF
697
698 ! Everything else, either via GauXC or manual XC computation
699 native_skala_evaluator = xc_section_uses_native_skala_grid(xc_section)
700 paw_composite_representation = .false.
701 IF (dft_control%use_gauxc .AND. (gapw .OR. gapw_xc)) THEN
702 paw_composite_representation = xc_section_uses_gauxc_model(xc_section) .AND. &
703 skala_gapw_representation(xc_section) == skala_gapw_paw_one_center .AND. &
704 gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
705 native_skala_evaluator = native_skala_evaluator .OR. &
706 paw_composite_representation
707 END IF
708 IF (dft_control%use_gauxc) THEN
709 IF (native_skala_evaluator) THEN
710 IF (gapw .AND. dft_control%low_spin_roks) THEN
711 CALL cp_abort(__location__, &
712 "GAPW low spin ROKS requires the molecular GauXC grid, not native Skala grids.")
713 END IF
714 CALL ensure_native_skala_grid_scope(xc_section)
715 gapw_representation = skala_gapw_representation(xc_section)
716 IF (gapw .OR. gapw_xc) THEN
717 gapw_all_all_electron = gauxc_gapw_all_all_electron(qs_kind_set)
718 gapw_all_pseudopotentials = gauxc_gapw_all_pseudopotentials(qs_kind_set)
719 gapw_mixed_core_representations = &
720 .NOT. gapw_all_all_electron .AND. .NOT. gapw_all_pseudopotentials
721 END IF
722 native_direct_valence_atom_composite = &
723 xc_section_uses_native_skala_grid(xc_section) .AND. &
724 native_skala_uses_atom_composite_grid(xc_section) .AND. &
725 ((.NOT. gapw .AND. .NOT. gapw_xc) .OR. &
726 (gapw_representation == skala_gapw_direct_valence .AND. &
727 gapw_all_pseudopotentials))
728 native_gapw_composite_reference = &
729 native_skala_gapw_composite_reference(xc_section) .AND. (gapw .OR. gapw_xc)
730 native_gapw_atom_composite_requested = &
731 native_skala_gapw_atom_composite_requested(xc_section) .AND. &
732 (gapw .OR. gapw_xc)
733 native_gapw_atom_composite_reference = &
734 native_gapw_atom_composite_requested .OR. paw_composite_representation
735 native_gapw_periodic_atom_composite = &
736 native_skala_uses_atom_composite_grid(xc_section)
737 IF (gapw_mixed_core_representations .AND. &
738 .NOT. native_gapw_periodic_atom_composite) THEN
739 CALL cp_abort(__location__, &
740 "Mixed all-electron and pseudopotential GAPW kinds require "// &
741 "NATIVE_GRID_LAYOUT ATOM_COMPOSITE so that their primitive "// &
742 "fields are combined before Skala feature construction.")
743 END IF
744 native_gapw_atom_composite_reference = &
745 native_gapw_atom_composite_reference .OR. &
746 (gapw_mixed_core_representations .AND. &
747 native_gapw_periodic_atom_composite)
748 native_gapw_all_electron_atom_composite = &
749 (gapw .OR. gapw_xc) .AND. &
750 gapw_all_all_electron .AND. &
751 (.NOT. native_gapw_composite_reference) .AND. &
752 native_gapw_periodic_atom_composite
753 native_gapw_atom_composite_reference = &
754 native_gapw_atom_composite_reference .OR. &
755 native_gapw_all_electron_atom_composite
756 IF (paw_composite_representation) THEN
757 IF (any(cell%perd /= 0)) THEN
758 native_gapw_composite_reference = &
759 .NOT. native_gapw_periodic_atom_composite
760 native_gapw_atom_composite_reference = &
761 native_gapw_periodic_atom_composite
762 ELSE
763 native_gapw_composite_reference = .false.
764 native_gapw_atom_composite_reference = .true.
765 END IF
766 END IF
767 IF (native_gapw_composite_reference .AND. native_gapw_atom_composite_reference) THEN
768 CALL cp_abort(__location__, &
769 "Select only one native-grid GAPW composite reference route.")
770 END IF
771 native_gapw_composite_direct_ao = native_gapw_composite_reference .AND. &
772 native_skala_gapw_composite_direct_ao(xc_section)
773 gauxc_section => get_gauxc_section(xc_section)
774 IF (ASSOCIATED(gauxc_section)) THEN
775 CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_DIAGNOSTICS", &
776 l_val=native_grid_diagnostics)
777 END IF
778 IF ((.NOT. do_kpoints) .AND. nimages /= 1) THEN
779 CALL cp_abort(__location__, &
780 "Native SKALA grid evaluation supports multiple images only "// &
781 "for k-point calculations.")
782 END IF
783 IF (do_kpoints) THEN
784 cpassert(ASSOCIATED(kpoints))
785 gauxc_section => get_gauxc_section(xc_section)
786 cpassert(ASSOCIATED(gauxc_section))
787 CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_USE_CUDA", l_val=native_grid_use_cuda)
788 IF (.NOT. native_grid_use_cuda) THEN
789 IF (para_env%mepos == 0 .AND. .NOT. native_grid_cpu_kpoints_warned) THEN
790 CALL cp_warn(__location__, &
791 "Native SKALA grid evaluation with k-points is using the CPU TorchScript "// &
792 "path. Use a mutually compatible BLAS, ScaLAPACK, OpenMP, and LibTorch "// &
793 "runtime stack, or select NATIVE_GRID_USE_CUDA T.")
794 native_grid_cpu_kpoints_warned = .true.
795 END IF
796 END IF
797 END IF
798 IF (dft_control%roks) THEN
799 cpabort("Native SKALA grid evaluation does not support ROKS.")
800 END IF
801 IF (dft_control%do_admm) THEN
802 cpabort("Native SKALA grid evaluation does not support ADMM.")
803 END IF
804 ! Force and stress rebuilds re-enter this path for derivatives and VXC only. The
805 ! autograd virial already contains the rebuilt XC volume term when qs_vxc_create
806 ! returns, so the externally reported energy must be restored to its SCF value.
807 native_skala_restore_exc = calculate_forces
808 IF (native_skala_restore_exc) THEN
809 native_skala_exc_scf = energy%exc
810 native_skala_total_scf = energy%total
811 END IF
812 IF (calculate_forces) THEN
813 ALLOCATE (native_skala_atom_force(3, natom))
814 native_skala_atom_force = 0.0_dp
815 IF (native_gapw_atom_composite_reference .OR. &
816 native_direct_valence_atom_composite) THEN
817 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
818 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
819 edisp=edisp, dispersion_env=qs_env%dispersion_env, &
820 just_energy=just_energy_xc, qs_env_external=qs_env, &
821 native_gapw_composite_override=native_gapw_composite_reference, &
822 native_skala_defer_to_atom_composite=.true.)
823 ELSE
824 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
825 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
826 edisp=edisp, dispersion_env=qs_env%dispersion_env, &
827 just_energy=just_energy_xc, &
828 native_skala_atom_force=native_skala_atom_force, &
829 qs_env_external=qs_env, &
830 native_gapw_composite_override=native_gapw_composite_reference)
831 END IF
832 IF (native_gapw_composite_reference .AND. .NOT. native_gapw_composite_direct_ao .AND. &
833 ASSOCIATED(v_rspace_new) .AND. ASSOCIATED(v_tau_rspace)) THEN
834 IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
835 output_unit = cp_logger_get_default_io_unit()
836 IF (output_unit > 0) THEN
837 DO iatom = 1, natom
838 WRITE (unit=output_unit, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
839 "SKALA_GPW| Native regular-grid atom force", iatom, &
840 native_skala_atom_force(:, iatom)
841 END DO
842 END IF
843 END IF
844 IF (use_virial) THEN
845 native_skala_composite_virial = 0.0_dp
846 CALL integrate_gapw_composite_vxc_gspace( &
847 qs_env, auxbas_pw_pool, v_rspace_new, v_tau_rspace, &
848 atom_force=native_skala_atom_force, &
849 strain_virial=native_skala_composite_virial, &
850 one_center_contraction=native_skala_composite_contraction, &
851 one_center_rho_contraction=native_skala_composite_rho_contraction, &
852 one_center_tau_contraction=native_skala_composite_tau_contraction)
853 IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
854 output_unit = cp_logger_get_default_io_unit()
855 IF (output_unit > 0) THEN
856 WRITE (unit=output_unit, fmt="(T2,A)") &
857 "SKALA_GPW| Composite-direction one-center virial"
858 DO iatom = 1, 3
859 WRITE (unit=output_unit, fmt="(T2,A,1X,3ES20.10)") &
860 "SKALA_GPW|", native_skala_composite_virial(iatom, :)
861 END DO
862 END IF
863 END IF
864 virial%pv_xc = virial%pv_xc + native_skala_composite_virial
865 ELSE
866 CALL integrate_gapw_composite_vxc_gspace( &
867 qs_env, auxbas_pw_pool, v_rspace_new, v_tau_rspace, &
868 atom_force=native_skala_atom_force, &
869 one_center_contraction=native_skala_composite_contraction, &
870 one_center_rho_contraction=native_skala_composite_rho_contraction, &
871 one_center_tau_contraction=native_skala_composite_tau_contraction)
872 END IF
873 END IF
874 IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
875 output_unit = cp_logger_get_default_io_unit()
876 IF (output_unit > 0) THEN
877 WRITE (unit=output_unit, fmt="(T2,A,1X,ES20.10)") &
878 "SKALA_GPW| Composite-direction one-center VXC contraction", &
879 native_skala_composite_contraction
880 WRITE (unit=output_unit, fmt="(T2,A,1X,ES20.10)") &
881 "SKALA_GPW| Composite-direction one-center rho contraction", &
882 native_skala_composite_rho_contraction
883 WRITE (unit=output_unit, fmt="(T2,A,1X,ES20.10)") &
884 "SKALA_GPW| Composite-direction one-center tau contraction", &
885 native_skala_composite_tau_contraction
886 DO iatom = 1, natom
887 WRITE (unit=output_unit, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
888 "SKALA_GPW| Native atom force", iatom, native_skala_atom_force(:, iatom)
889 END DO
890 END IF
891 END IF
892 cpassert(ASSOCIATED(force))
893 cpassert(ASSOCIATED(atomic_kind_set))
894 CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
895 DO iatom = 1, natom
896 ikind = kind_of(iatom)
897 atom_a = atom_of_kind(iatom)
898 force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + &
899 native_skala_atom_force(:, iatom)
900 END DO
901 DEALLOCATE (atom_of_kind, kind_of, native_skala_atom_force)
902 ELSE
903 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
904 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
905 edisp=edisp, dispersion_env=qs_env%dispersion_env, &
906 just_energy=just_energy_xc, qs_env_external=qs_env, &
907 native_gapw_composite_override=native_gapw_composite_reference, &
908 native_skala_defer_to_atom_composite= &
909 (native_gapw_atom_composite_reference .OR. &
910 native_direct_valence_atom_composite))
911 IF (native_gapw_composite_reference .AND. .NOT. native_gapw_composite_direct_ao .AND. &
912 ASSOCIATED(v_rspace_new) .AND. ASSOCIATED(v_tau_rspace)) THEN
913 CALL integrate_gapw_composite_vxc_gspace( &
914 qs_env, auxbas_pw_pool, v_rspace_new, v_tau_rspace, &
915 one_center_contraction=native_skala_composite_contraction, &
916 one_center_rho_contraction=native_skala_composite_rho_contraction, &
917 one_center_tau_contraction=native_skala_composite_tau_contraction)
918 IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
919 output_unit = cp_logger_get_default_io_unit()
920 IF (output_unit > 0) THEN
921 WRITE (unit=output_unit, fmt="(T2,A,1X,ES20.10)") &
922 "SKALA_GPW| Composite-direction one-center VXC contraction", &
923 native_skala_composite_contraction
924 WRITE (unit=output_unit, fmt="(T2,A,1X,ES20.10)") &
925 "SKALA_GPW| Composite-direction one-center rho contraction", &
926 native_skala_composite_rho_contraction
927 WRITE (unit=output_unit, fmt="(T2,A,1X,ES20.10)") &
928 "SKALA_GPW| Composite-direction one-center tau contraction", &
929 native_skala_composite_tau_contraction
930 END IF
931 END IF
932 END IF
933 END IF
934 IF (native_skala_restore_exc) energy%exc = native_skala_exc_scf
935 IF (native_direct_valence_atom_composite) THEN
936 cpassert(ASSOCIATED(v_rspace_new))
937 cpassert(ASSOCIATED(v_tau_rspace))
938 CALL calculate_vxc_atom( &
939 qs_env, just_energy_xc, energy%exc1, xc_section_external=xc_section, &
940 calculate_forces=calculate_forces, composite_vxc_rho=v_rspace_new, &
941 composite_vxc_tau=v_tau_rspace, direct_valence_atom_grid=.true.)
942 energy%exc = 0.0_dp
943 ELSE IF ((gapw .OR. gapw_xc) .AND. .NOT. native_gapw_composite_reference) THEN
944 IF (native_gapw_atom_composite_reference) THEN
945 cpassert(ASSOCIATED(v_rspace_new))
946 cpassert(ASSOCIATED(v_tau_rspace))
947 CALL calculate_vxc_atom( &
948 qs_env, just_energy_xc, energy%exc1, xc_section_external=xc_section, &
949 calculate_forces=calculate_forces, composite_vxc_rho=v_rspace_new, &
950 composite_vxc_tau=v_tau_rspace, atom_composite_grid=.true.)
951 energy%exc = 0.0_dp
952 ELSE
953 CALL calculate_vxc_atom(qs_env, just_energy_xc, energy%exc1, &
954 xc_section_external=xc_section, &
955 calculate_forces=calculate_forces)
956 END IF
957 END IF
958 IF (edisp /= 0.0_dp) energy%dispersion = edisp
959 IF (qs_env%requires_matrix_vxc .AND. ASSOCIATED(v_rspace_new)) THEN
960 IF (do_kpoints) THEN
961 CALL compute_matrix_vxc_kp(qs_env=qs_env, v_rspace=v_rspace_new, &
962 matrix_vxc_kp=matrix_vxc_kp, &
963 gapw_full_basis=native_gapw_composite_direct_ao)
964 CALL set_ks_env(ks_env, matrix_vxc_kp=matrix_vxc_kp)
965 ELSE
966 CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace_new, &
967 matrix_vxc=matrix_vxc, &
968 gapw_full_basis=native_gapw_composite_direct_ao)
969 CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
970 END IF
971 END IF
972 ELSE
973 use_gauxc_matrix = .true.
974 CALL apply_gauxc(qs_env, xc_section, calculate_forces)
975 IF (gapw_xc .OR. (gapw .AND. gauxc_gapw_has_paw_pseudopotentials(qs_kind_set))) THEN
976 CALL calculate_vxc_atom(qs_env, just_energy_xc, energy%exc1, &
977 xc_section_external=xc_section, &
978 calculate_forces=calculate_forces)
979 END IF
980 END IF
981 ELSE
982 CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
983 vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
984 edisp=edisp, dispersion_env=qs_env%dispersion_env, &
985 just_energy=just_energy_xc)
986 IF (edisp /= 0.0_dp) energy%dispersion = edisp
987 IF (qs_env%requires_matrix_vxc .AND. ASSOCIATED(v_rspace_new)) THEN
988 CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace_new, matrix_vxc=matrix_vxc)
989 CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
990 END IF
991
992 IF (gapw .OR. gapw_xc) THEN
993 CALL calculate_vxc_atom(qs_env, just_energy_xc, energy%exc1, &
994 xc_section_external=xc_section, &
995 calculate_forces=calculate_forces)
996 END IF
997 END IF
998 END IF
999
1000 IF (gapw .AND. dft_control%qs_control%cdft) THEN
1001 ALLOCATE (cdft_gapw_values(SIZE(cdft_control%group)))
1002 IF (cdft_control%atomic_charges) THEN
1003 ALLOCATE (cdft_charge_correction(natom, nspins), &
1004 cdft_electronic_charge(cdft_control%natoms, nspins))
1005 CALL gapw_cdft_one_center(qs_env, just_energy, calculate_forces, &
1006 cdft_gapw_values, cdft_charge_correction)
1007 DO ispin = 1, nspins
1008 DO iatom = 1, cdft_control%natoms
1009 cdft_electronic_charge(iatom, ispin) = &
1010 pw_integral_ab(cdft_control%charge(iatom), rho_r(ispin), local_only=.true.)
1011 END DO
1012 END DO
1013 CALL para_env%sum(cdft_electronic_charge)
1014 DO ispin = 1, nspins
1015 DO iatom = 1, cdft_control%natoms
1016 cdft_electronic_charge(iatom, ispin) = &
1017 cdft_electronic_charge(iatom, ispin) + &
1018 cdft_charge_correction(cdft_control%atoms(iatom), ispin)
1019 END DO
1020 END DO
1021 ELSE
1022 CALL gapw_cdft_one_center(qs_env, just_energy, calculate_forces, cdft_gapw_values)
1023 END IF
1024 cdft_control%value = cdft_control%value + cdft_gapw_values
1025 energy%cdft = energy%cdft + dot_product(cdft_control%strength, cdft_gapw_values)
1026 IF (cdft_control%atomic_charges) THEN
1027 CALL cdft_constraint_print(qs_env, cdft_electronic_charge)
1028 DEALLOCATE (cdft_charge_correction, cdft_electronic_charge)
1029 ELSE
1030 CALL cdft_constraint_print(qs_env, cdft_electronic_charge)
1031 END IF
1032 DEALLOCATE (cdft_gapw_values)
1033 END IF
1034
1035 IF (dft_control%mtlr_dft_with_perturbation .AND. .NOT. just_energy) THEN
1036 IF (do_adiabatic_rescaling .OR. use_gauxc_matrix) THEN
1037 cpabort("MTLR requires an XC potential on the real-space grid.")
1038 END IF
1039 IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
1040 ALLOCATE (matrix_vhxc(nspins))
1041 DO ispin = 1, nspins
1042 NULLIFY (matrix_vhxc(ispin)%matrix)
1043 ALLOCATE (matrix_vhxc(ispin)%matrix)
1044 CALL dbcsr_create(matrix_vhxc(ispin)%matrix, &
1045 template=ks_matrix(1, 1)%matrix)
1046 CALL cp_dbcsr_alloc_block_from_nbl(matrix_vhxc(ispin)%matrix, sab_orb)
1047 CALL dbcsr_set(matrix_vhxc(ispin)%matrix, 0.0_dp)
1048 END DO
1049 CALL set_ks_env(ks_env, matrix_vhxc=matrix_vhxc)
1050 ELSE
1051 cpassert(SIZE(matrix_vhxc) == nspins)
1052 DO ispin = 1, nspins
1053 cpassert(ASSOCIATED(matrix_vhxc(ispin)%matrix))
1054 END DO
1055 END IF
1056 DO ispin = 1, nspins
1057 CALL auxbas_pw_pool%create_pw(v_hxc_rspace(ispin))
1058 CALL pw_zero(v_hxc_rspace(ispin))
1059 CALL pw_axpy(v_hartree_rspace, v_hxc_rspace(ispin))
1060 IF (ASSOCIATED(v_rspace_new)) THEN
1061 CALL pw_axpy(v_rspace_new(ispin), &
1062 v_hxc_rspace(ispin), &
1063 v_rspace_new(ispin)%pw_grid%dvol)
1064 END IF
1065 CALL dbcsr_set(matrix_vhxc(ispin)%matrix, 0.0_dp)
1066 CALL integrate_v_rspace(v_hxc_rspace(ispin), &
1067 hmat=matrix_vhxc(ispin), &
1068 qs_env=qs_env, &
1069 calculate_forces=.false.)
1070 CALL auxbas_pw_pool%give_back_pw(v_hxc_rspace(ispin))
1071 END DO
1072 END IF
1073
1074 ! set hartree and xc potentials for use in Harris method
1075 IF (qs_env%harris_method) THEN
1076 CALL get_qs_env(qs_env, harris_env=harris_env)
1077 CALL harris_set_potentials(harris_env, v_hartree_rspace, v_rspace_new)
1078 END IF
1079
1080 NULLIFY (rho_struct)
1081 IF (use_virial .AND. calculate_forces) THEN
1082 virial%pv_exc = virial%pv_exc - virial%pv_xc
1083 virial%pv_virial = virial%pv_virial - virial%pv_xc
1084 END IF
1085
1086 ! *** Add Hartree-Fock contribution if required ***
1087 hfx_sections => section_vals_get_subs_vals(xc_section, "HF")
1088 CALL section_vals_get(hfx_sections, explicit=do_hfx)
1089
1090 ace_active = .false.
1091 ace_rebuild_frequency = 1
1092
1093 IF (do_hfx) THEN
1094 ace_section => section_vals_get_subs_vals(hfx_sections, "ACE")
1095 IF (ASSOCIATED(ace_section)) THEN
1096 CALL section_vals_val_get(ace_section, "ACTIVE", l_val=ace_active)
1097 CALL section_vals_val_get(ace_section, "REBUILD_FREQUENCY", i_val=ace_rebuild_frequency)
1098 END IF
1099 END IF
1100
1101 IF (do_hfx) THEN
1102 IF (dokp) THEN
1103 IF (ace_active) THEN
1104 cpabort("ACE-HFX for k-points is not implemented yet")
1105 ELSE
1106 CALL hfx_ks_matrix_kp(qs_env, ks_matrix, energy, calculate_forces)
1107 END IF
1108
1109 ELSE
1110 ! ext_xc_section may contain a hfx section
1111 IF (ace_active) THEN
1112 CALL hfx_ace_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, &
1113 just_energy, v_rspace_new, v_tau_rspace, &
1114 ace_rebuild_frequency, ext_xc_section=xc_section)
1115 ELSE
1116 CALL hfx_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, &
1117 just_energy, v_rspace_new, v_tau_rspace, ext_xc_section=xc_section)
1118 END IF
1119 END IF
1120 END IF !do_hfx
1121
1122 IF (do_ppl .AND. calculate_forces) THEN
1123 cpassert(.NOT. gapw)
1124 DO ispin = 1, nspins
1125 CALL integrate_ppl_rspace(rho_r(ispin), qs_env)
1126 END DO
1127 END IF
1128
1129 ! The atom-composite evaluator differentiates its atom-centered NLCC field directly.
1130 IF (ASSOCIATED(rho_nlcc) .AND. calculate_forces .AND. &
1131 .NOT. native_gapw_atom_composite_reference .AND. &
1132 .NOT. native_direct_valence_atom_composite) THEN
1133 DO ispin = 1, nspins
1134 CALL integrate_rho_nlcc(v_rspace_new(ispin), qs_env)
1135 IF (dft_control%do_admm) CALL integrate_rho_nlcc(v_rspace_new_aux_fit(ispin), qs_env)
1136 END DO
1137 END IF
1138
1139 ! calculate KG correction
1140 IF (dft_control%qs_control%do_kg .AND. just_energy) THEN
1141
1142 cpassert(nimages == 1)
1143 ksmat => ks_matrix(:, 1)
1144 CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces, do_kernel=.false.)
1145
1146 ! subtract kg corr from the total energy
1147 energy%exc = energy%exc - ekin_mol
1148
1149 END IF
1150
1151 ! *** Single atom contributions ***
1152 IF (.NOT. just_energy) THEN
1153 IF (calculate_forces) THEN
1154 ! Getting nuclear force contribution from the core charge density
1155 IF ((poisson_env%parameters%solver == pw_poisson_implicit) .AND. &
1156 (poisson_env%parameters%dielectric_params%dielec_core_correction)) THEN
1157 block
1158 TYPE(pw_r3d_rs_type) :: v_minus_veps
1159 CALL auxbas_pw_pool%create_pw(v_minus_veps)
1160 CALL pw_copy(v_hartree_rspace, v_minus_veps)
1161 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_minus_veps, -v_hartree_rspace%pw_grid%dvol)
1162 CALL integrate_v_core_rspace(v_minus_veps, qs_env)
1163 CALL auxbas_pw_pool%give_back_pw(v_minus_veps)
1164 END block
1165 ELSE
1166 CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
1167 END IF
1168 ! SCCS with core fill: the dielectric function depends on rho_elec + rho_fill(R),
1169 ! which yields an additional nuclear force contribution from v_sccs. The fill
1170 ! Gaussians enter integrate_v_core_rspace with the sign convention of rho_core
1171 ! (-ccore*exp(-calpha*r**2)), hence the negated potential
1172 IF (dft_control%do_sccs) THEN
1173 IF (dft_control%sccs_control%core_fill) THEN
1174 block
1175 REAL(kind=dp), DIMENSION(:), ALLOCATABLE :: calpha, ccore
1176 TYPE(pw_r3d_rs_type) :: v_minus_vsccs
1177 CALL sccs_core_fill_coefficients(qs_env, calpha, ccore)
1178 CALL auxbas_pw_pool%create_pw(v_minus_vsccs)
1179 CALL pw_copy(v_sccs_rspace, v_minus_vsccs)
1180 CALL pw_scale(v_minus_vsccs, -1.0_dp)
1181 CALL integrate_v_core_rspace(v_minus_vsccs, qs_env, calpha=calpha, ccore=ccore)
1182 CALL auxbas_pw_pool%give_back_pw(v_minus_vsccs)
1183 DEALLOCATE (calpha, ccore)
1184 END block
1185 END IF
1186 END IF
1187 END IF
1188
1189 IF (.NOT. do_hfx) THEN
1190 ! Initialize the Kohn-Sham matrix with the core Hamiltonian matrix
1191 ! (sets ks sparsity equal to matrix_h sparsity)
1192 DO ispin = 1, nspins
1193 DO img = 1, nimages
1194 CALL dbcsr_get_info(ks_matrix(ispin, img)%matrix, name=name) ! keep the name
1195 CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix, name=name)
1196 END DO
1197 END DO
1198 ! imaginary part if required
1199 IF (qs_env%run_rtp) THEN
1200 IF (dft_control%rtp_control%velocity_gauge) THEN
1201 cpassert(ASSOCIATED(matrix_h_im))
1202 cpassert(ASSOCIATED(ks_matrix_im))
1203 DO ispin = 1, nspins
1204 DO img = 1, nimages
1205 CALL dbcsr_get_info(ks_matrix_im(ispin, img)%matrix, name=name) ! keep the name
1206 CALL dbcsr_copy(ks_matrix_im(ispin, img)%matrix, matrix_h_im(1, img)%matrix, name=name)
1207 END DO
1208 END DO
1209 END IF
1210 END IF
1211 END IF
1212
1213 IF (use_virial .AND. calculate_forces) THEN
1214 pv_loc = virial%pv_virial
1215 END IF
1216 ! sum up potentials and integrate
1217 ! Pointing my_rho to the density matrix rho_ao
1218 my_rho => rho_ao
1219
1220 CALL sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, &
1221 v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, &
1222 v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, &
1223 cdft_control, calculate_forces)
1224
1225 IF (use_gauxc_matrix) THEN
1226 IF (dokp) THEN
1227 CALL get_qs_env(qs_env=qs_env, matrix_vxc_kp=matrix_vxc_kp)
1228 cpassert(ASSOCIATED(matrix_vxc_kp))
1229 DO ispin = 1, nspins
1230 DO img = 1, nimages
1231 CALL dbcsr_add(ks_matrix(ispin, img)%matrix, matrix_vxc_kp(ispin, img)%matrix, &
1232 1.0_dp, 1.0_dp)
1233 END DO
1234 END DO
1235 ELSE
1236 CALL get_qs_env(qs_env=qs_env, matrix_vxc=matrix_vxc)
1237 cpassert(ASSOCIATED(matrix_vxc))
1238 cpassert(nimages == 1)
1239 DO ispin = 1, nspins
1240 CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, matrix_vxc(ispin)%matrix, 1.0_dp, 1.0_dp)
1241 END DO
1242 END IF
1243 END IF
1244
1245 IF (gapw .OR. gapw_xc) THEN
1246 IF (calculate_forces) THEN
1247 IF (gapw_xc) THEN
1248 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
1249 ELSE
1250 CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
1251 END IF
1252 NULLIFY (rho1)
1253 IF (dft_control%use_gauxc .AND. (gapw .OR. gapw_xc) .AND. &
1254 .NOT. native_skala_evaluator) THEN
1255 gauxc_model_none = .false.
1256 gauxc_section => get_gauxc_section(xc_section)
1257 IF (ASSOCIATED(gauxc_section)) THEN
1258 CALL section_vals_val_get(gauxc_section, "MODEL", c_val=gauxc_model_name)
1259 gauxc_model_name = adjustl(gauxc_model_name)
1260 CALL uppercase(gauxc_model_name)
1261 gauxc_model_none = (trim(gauxc_model_name) == "" .OR. &
1262 trim(gauxc_model_name) == "NONE")
1263 END IF
1264 IF (gauxc_model_none .AND. &
1265 (gapw_xc .OR. gauxc_gapw_has_paw_pseudopotentials(qs_kind_set))) THEN
1266 CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section)
1267 END IF
1268 ELSE
1269 CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section)
1270 END IF
1271 !
1272 IF (dft_control%do_admm) THEN
1273 CALL get_qs_env(qs_env, admm_env=admm_env)
1274 xc_section => admm_env%xc_section_aux
1275 CALL get_admm_env(admm_env, rho_aux_fit=rho_struct)
1276 vscale = 1.0_dp
1277 IF (admm_env%do_admmp) THEN
1278 vscale = admm_env%gsi(1)**2
1279 ELSE IF (admm_env%do_admms) THEN
1280 vscale = admm_env%gsi(1)**(2.0_dp/3.0_dp)
1281 END IF
1282 CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section, force_scale=vscale)
1283 END IF
1284 END IF
1285 END IF
1286
1287 IF (use_virial .AND. calculate_forces) THEN
1288 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1289 END IF
1290 IF (dft_control%qs_control%do_kg) THEN
1291 cpassert(nimages == 1)
1292 ksmat => ks_matrix(:, 1)
1293
1294 IF (use_virial .AND. calculate_forces) THEN
1295 pv_loc = virial%pv_virial
1296 END IF
1297
1298 CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces, do_kernel=.false.)
1299 ! subtract kg corr from the total energy
1300 energy%exc = energy%exc - ekin_mol
1301
1302 ! virial corrections
1303 IF (use_virial .AND. calculate_forces) THEN
1304
1305 ! Integral contribution
1306 virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1307
1308 ! GGA contribution
1309 virial%pv_exc = virial%pv_exc + virial%pv_xc
1310 virial%pv_virial = virial%pv_virial + virial%pv_xc
1311 virial%pv_xc = 0.0_dp
1312 END IF
1313 END IF
1314
1315 ELSE
1316 ! Atom-composite XC also creates potential buffers for energy-only calls.
1317 ! These calls skip sum_up_and_integrate, which normally returns the buffers.
1318 IF (ASSOCIATED(v_rspace_new)) THEN
1319 DO ispin = 1, SIZE(v_rspace_new)
1320 CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
1321 END DO
1322 DEALLOCATE (v_rspace_new)
1323 END IF
1324 IF (ASSOCIATED(v_tau_rspace)) THEN
1325 DO ispin = 1, SIZE(v_tau_rspace)
1326 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1327 END DO
1328 DEALLOCATE (v_tau_rspace)
1329 END IF
1330 IF (do_hfx) THEN
1331 IF (.false.) THEN
1332 cpwarn("KS matrix no longer correct. Check possible problems with property calculations!")
1333 END IF
1334 END IF
1335 END IF ! .NOT. just energy
1336 IF (dft_control%qs_control%ddapc_explicit_potential) THEN
1337 CALL auxbas_pw_pool%give_back_pw(v_spin_ddapc_rest_r)
1338 DEALLOCATE (v_spin_ddapc_rest_r)
1339 END IF
1340
1341 IF (calculate_forces .AND. dft_control%qs_control%cdft) THEN
1342 IF (.NOT. cdft_control%transfer_pot) THEN
1343 DO iatom = 1, SIZE(cdft_control%group)
1344 CALL auxbas_pw_pool%give_back_pw(cdft_control%group(iatom)%weight)
1345 DEALLOCATE (cdft_control%group(iatom)%weight)
1346 END DO
1347 IF (cdft_control%atomic_charges) THEN
1348 DO iatom = 1, cdft_control%natoms
1349 CALL auxbas_pw_pool%give_back_pw(cdft_control%charge(iatom))
1350 END DO
1351 DEALLOCATE (cdft_control%charge)
1352 END IF
1353 IF (cdft_control%type == outer_scf_becke_constraint .AND. &
1354 cdft_control%becke_control%cavity_confine) THEN
1355 IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN
1356 CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
1357 ELSE
1358 DEALLOCATE (cdft_control%becke_control%cavity_mat)
1359 END IF
1360 ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1361 IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) THEN
1362 CALL auxbas_pw_pool%give_back_pw(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)
1363 END IF
1364 END IF
1365 IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment)
1366 cdft_control%save_pot = .false.
1367 cdft_control%need_pot = .true.
1368 cdft_control%external_control = .false.
1369 END IF
1370 END IF
1371
1372 IF (dft_control%do_sccs) THEN
1373 CALL auxbas_pw_pool%give_back_pw(v_sccs_rspace)
1374 DEALLOCATE (v_sccs_rspace)
1375 END IF
1376
1377 IF (gapw) THEN
1378 IF (dft_control%apply_external_potential) THEN
1379 ! Integrals of the Hartree potential with g0_soft
1380 CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, &
1381 v_qmmm=vee, scale=-1.0_dp)
1382 END IF
1383 CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces)
1384 ! Place Vh_1c_gg_integrals after integrate_vhg0_rspace for CNEO calculations
1385 ! because vhg0 integral is needed to build the complete nuclear equation
1386 CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c, local_rho_set=local_rho_set)
1387 CALL vh_1c_gg_integrals(qs_env, energy%hartree_1c, ecoul_1c, local_rho_set, para_env, tddft=.false., &
1388 core_2nd=.false.)
1389 ! CNEO quantum nuclear core energy (kinetic + Z*erfc(r)/r potential from classical nuclei)
1390 energy%core_cneo = 0.0_dp
1391 IF (ASSOCIATED(local_rho_set%rhoz_cneo_set)) THEN
1392 DO iatom = 1, SIZE(local_rho_set%rhoz_cneo_set)
1393 energy%core_cneo = energy%core_cneo + local_rho_set%rhoz_cneo_set(iatom)%e_core
1394 END DO
1395 END IF
1396 END IF
1397
1398 IF (gapw .OR. gapw_xc) THEN
1399 ! Single atom contributions in the KS matrix ***
1400 CALL update_ks_atom(qs_env, ks_matrix, rho_ao, calculate_forces)
1401 IF (dft_control%do_admm) THEN
1402 !Single atom contribution to the AUX matrices
1403 !Note: also update ks_aux_fit matrix in case of rtp
1404 CALL admm_update_ks_atom(qs_env, calculate_forces)
1405 END IF
1406 END IF
1407
1408 !Calculation of Mulliken restraint, if requested
1409 CALL qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, &
1410 ks_matrix, matrix_s, rho, mulliken_order_p)
1411
1412 ! Add DFT+U contribution, if requested
1413 IF (dft_control%dft_plus_u) THEN
1414 IF (just_energy) THEN
1415 CALL plus_u(qs_env=qs_env)
1416 ELSE
1417 CALL plus_u(qs_env=qs_env, matrix_h=ks_matrix)
1418 END IF
1419 ELSE
1420 energy%dft_plus_u = 0.0_dp
1421 END IF
1422
1423 ! At this point the ks matrix should be up to date, filter it if requested
1424 DO ispin = 1, nspins
1425 DO img = 1, nimages
1426 CALL dbcsr_filter(ks_matrix(ispin, img)%matrix, &
1427 dft_control%qs_control%eps_filter_matrix)
1428 END DO
1429 END DO
1430
1431 !** merge the auxiliary KS matrix and the primary one
1432 IF (dft_control%do_admm_mo) THEN
1433 IF (qs_env%run_rtp) THEN
1434 CALL rtp_admm_merge_ks_matrix(qs_env)
1435 ELSE
1436 CALL admm_mo_merge_ks_matrix(qs_env)
1437 END IF
1438 ELSE IF (dft_control%do_admm_dm) THEN
1439 CALL admm_dm_merge_ks_matrix(qs_env)
1440 END IF
1441
1442 ! External field (nonperiodic case)
1443 CALL qs_efield_local_operator(qs_env, just_energy, calculate_forces)
1444
1445 ! Right now we can compute the orbital derivative here, as it depends currently only on the available
1446 ! Kohn-Sham matrix. This might change in the future, in which case more pieces might need to be assembled
1447 ! from this routine, notice that this part of the calculation in not linear scaling
1448 ! right now this operation is only non-trivial because of occupation numbers and the restricted keyword
1449 IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy .AND. .NOT. qs_env%run_rtp) THEN
1450 CALL get_qs_env(qs_env, mo_derivs=mo_derivs)
1451 cpassert(nimages == 1)
1452 ksmat => ks_matrix(:, 1)
1453 CALL calc_mo_derivatives(qs_env, ksmat, mo_derivs)
1454 END IF
1455
1456 ! ADMM overlap forces
1457 IF (calculate_forces .AND. dft_control%do_admm) THEN
1458 IF (dokp) THEN
1459 CALL calc_admm_ovlp_forces_kp(qs_env)
1460 ELSE
1461 CALL calc_admm_ovlp_forces(qs_env)
1462 END IF
1463 END IF
1464
1465 ! deal with low spin roks
1466 CALL low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
1467 calculate_forces, auxbas_pw_pool)
1468
1469 ! deal with sic on explicit orbitals
1470 CALL sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, &
1471 calculate_forces, auxbas_pw_pool)
1472
1473 ! Periodic external field
1474 CALL qs_efield_berry_phase(qs_env, just_energy, calculate_forces)
1475
1476 ! adds s2_restraint energy and orbital derivatives
1477 CALL qs_ks_s2_restraint(dft_control, qs_env, matrix_s, &
1478 energy, calculate_forces, just_energy)
1479
1480 IF (do_ppl) THEN
1481 ! update core energy for grid based local pseudopotential
1482 ecore_ppl = 0._dp
1483 DO ispin = 1, nspins
1484 ecore_ppl = ecore_ppl + pw_integral_ab(vppl_rspace, rho_r(ispin))
1485 END DO
1486 energy%core = energy%core + ecore_ppl
1487 END IF
1488
1489 IF (lrigpw) THEN
1490 ! update core energy for ppl_ri method
1491 CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density)
1492 IF (lri_env%ppl_ri) THEN
1493 ecore_ppl = 0._dp
1494 DO ispin = 1, nspins
1495 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1496 CALL v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl)
1497 END DO
1498 energy%core = energy%core + ecore_ppl
1499 END IF
1500 END IF
1501
1502 ! Sum all energy terms to obtain the total energy
1503 IF (native_skala_restore_exc) THEN
1504 energy%total = native_skala_total_scf
1505 IF (abnormal_value(energy%total)) THEN
1506 CALL cp_abort(__location__, &
1507 "SCF energy restored by native Skala is an abnormal value (NaN/Inf).")
1508 END IF
1509 ELSE
1510 CALL check_sum_energies(energy)
1511 END IF
1512
1513 ! Print detailed energy
1514 IF (my_print) THEN
1515 CALL print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
1516 END IF
1517
1518 CALL timestop(handle)
1519
1520 END SUBROUTINE qs_ks_build_kohn_sham_matrix
1521
1522! **************************************************************************************************
1523!> \brief ...
1524!> \param rho_tot_gspace ...
1525!> \param qs_env ...
1526!> \param rho ...
1527!> \param skip_nuclear_density ...
1528! **************************************************************************************************
1529 SUBROUTINE calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
1530 TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_tot_gspace
1531 TYPE(qs_environment_type), POINTER :: qs_env
1532 TYPE(qs_rho_type), POINTER :: rho
1533 LOGICAL, INTENT(IN), OPTIONAL :: skip_nuclear_density
1534
1535 INTEGER :: ispin
1536 LOGICAL :: my_skip
1537 TYPE(dft_control_type), POINTER :: dft_control
1538 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1539 TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
1540 TYPE(qs_charges_type), POINTER :: qs_charges
1541
1542 my_skip = .false.
1543 IF (PRESENT(skip_nuclear_density)) my_skip = skip_nuclear_density
1544
1545 CALL qs_rho_get(rho, rho_g=rho_g)
1546 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
1547
1548 IF (.NOT. my_skip) THEN
1549 NULLIFY (rho_core)
1550 CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1551 IF (dft_control%qs_control%gapw) THEN
1552 NULLIFY (rho0_s_gs, rhoz_cneo_s_gs)
1553 CALL get_qs_env(qs_env=qs_env, rho0_s_gs=rho0_s_gs, rhoz_cneo_s_gs=rhoz_cneo_s_gs)
1554 cpassert(ASSOCIATED(rho0_s_gs))
1555 CALL pw_copy(rho0_s_gs, rho_tot_gspace)
1556 IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
1557 CALL pw_axpy(rhoz_cneo_s_gs, rho_tot_gspace)
1558 END IF
1559 IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
1560 CALL pw_axpy(rho_core, rho_tot_gspace)
1561 END IF
1562 ELSE
1563 CALL pw_copy(rho_core, rho_tot_gspace)
1564 END IF
1565 DO ispin = 1, dft_control%nspins
1566 CALL pw_axpy(rho_g(ispin), rho_tot_gspace)
1567 END DO
1568 CALL get_qs_env(qs_env=qs_env, qs_charges=qs_charges)
1569 qs_charges%total_rho_gspace = pw_integrate_function(rho_tot_gspace, isign=-1)
1570 ELSE
1571 DO ispin = 1, dft_control%nspins
1572 CALL pw_axpy(rho_g(ispin), rho_tot_gspace)
1573 END DO
1574 END IF
1575
1576 END SUBROUTINE calc_rho_tot_gspace
1577
1578! **************************************************************************************************
1579!> \brief compute MO derivatives
1580!> \param qs_env the qs_env to update
1581!> \param ks_matrix ...
1582!> \param mo_derivs ...
1583!> \par History
1584!> 01.2014 created, transferred from qs_ks_build_kohn_sham_matrix in
1585!> separate subroutine
1586!> \author Dorothea Golze
1587! **************************************************************************************************
1588 SUBROUTINE calc_mo_derivatives(qs_env, ks_matrix, mo_derivs)
1589 TYPE(qs_environment_type), POINTER :: qs_env
1590 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, mo_derivs
1591
1592 INTEGER :: ispin
1593 LOGICAL :: uniform_occupation
1594 REAL(kind=dp), DIMENSION(:), POINTER :: occupation_numbers
1595 TYPE(cp_fm_type), POINTER :: mo_coeff
1596 TYPE(dbcsr_type) :: mo_derivs2_tmp1, mo_derivs2_tmp2
1597 TYPE(dbcsr_type), POINTER :: mo_coeff_b
1598 TYPE(dft_control_type), POINTER :: dft_control
1599 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
1600
1601 NULLIFY (dft_control, mo_array, mo_coeff, mo_coeff_b, occupation_numbers)
1602
1603 CALL get_qs_env(qs_env, &
1604 dft_control=dft_control, &
1605 mos=mo_array)
1606
1607 DO ispin = 1, SIZE(mo_derivs)
1608
1609 CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, &
1610 mo_coeff_b=mo_coeff_b, occupation_numbers=occupation_numbers)
1611 CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin)%matrix, mo_coeff_b, &
1612 0.0_dp, mo_derivs(ispin)%matrix)
1613
1614 IF (dft_control%restricted) THEN
1615 ! only the first mo_set are actual variables, but we still need both
1616 cpassert(ispin == 1)
1617 cpassert(SIZE(mo_array) == 2)
1618 ! use a temporary array with the same size as the first spin for the second spin
1619
1620 ! uniform_occupation is needed for this case, otherwise we can not
1621 ! reconstruct things in ot, since we irreversibly sum
1622 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
1623 cpassert(uniform_occupation)
1624 CALL get_mo_set(mo_set=mo_array(2), uniform_occupation=uniform_occupation)
1625 cpassert(uniform_occupation)
1626
1627 ! The beta-spin might have fewer orbitals than alpa-spin...
1628 ! create temporary matrices with beta_nmo columns
1629 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff_b)
1630 CALL dbcsr_create(mo_derivs2_tmp1, template=mo_coeff_b)
1631
1632 ! calculate beta derivatives
1633 CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(2)%matrix, mo_coeff_b, 0.0_dp, mo_derivs2_tmp1)
1634
1635 ! create larger matrix with alpha_nmo columns
1636 CALL dbcsr_create(mo_derivs2_tmp2, template=mo_derivs(1)%matrix)
1637 CALL dbcsr_set(mo_derivs2_tmp2, 0.0_dp)
1638
1639 ! copy into larger matrix, fills the first beta_nmo columns
1640 CALL dbcsr_copy_columns_hack(mo_derivs2_tmp2, mo_derivs2_tmp1, &
1641 mo_array(2)%nmo, 1, 1, &
1642 para_env=mo_array(1)%mo_coeff%matrix_struct%para_env, &
1643 blacs_env=mo_array(1)%mo_coeff%matrix_struct%context)
1644
1645 ! add beta contribution to alpa mo_derivs
1646 CALL dbcsr_add(mo_derivs(1)%matrix, mo_derivs2_tmp2, 1.0_dp, 1.0_dp)
1647 CALL dbcsr_release(mo_derivs2_tmp1)
1648 CALL dbcsr_release(mo_derivs2_tmp2)
1649 END IF
1650 END DO
1651
1652 IF (dft_control%do_admm_mo) THEN
1653 CALL calc_admm_mo_derivatives(qs_env, mo_derivs)
1654 END IF
1655
1656 END SUBROUTINE calc_mo_derivatives
1657
1658! **************************************************************************************************
1659!> \brief updates the Kohn Sham matrix of the given qs_env (facility method)
1660!> \param qs_env the qs_env to update
1661!> \param calculate_forces if true calculate the quantities needed
1662!> to calculate the forces. Defaults to false.
1663!> \param just_energy if true updates the energies but not the
1664!> ks matrix. Defaults to false
1665!> \param print_active ...
1666!> \par History
1667!> 4.2002 created [fawzi]
1668!> 8.2014 kpoints [JGH]
1669!> 10.2014 refractored [Ole Schuett]
1670!> \author Fawzi Mohamed
1671! **************************************************************************************************
1672 SUBROUTINE qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, &
1673 print_active)
1674 TYPE(qs_environment_type), POINTER :: qs_env
1675 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces, just_energy, &
1676 print_active
1677
1678 CHARACTER(LEN=*), PARAMETER :: routinen = 'qs_ks_update_qs_env'
1679
1680 INTEGER :: handle, unit_nr
1681 LOGICAL :: c_forces, do_rebuild, energy_only, &
1682 forces_up_to_date, potential_changed, &
1683 rho_changed, s_mstruct_changed
1684 TYPE(qs_ks_env_type), POINTER :: ks_env
1685
1686 NULLIFY (ks_env)
1687 unit_nr = cp_logger_get_default_io_unit()
1688
1689 c_forces = .false.
1690 energy_only = .false.
1691 IF (PRESENT(just_energy)) energy_only = just_energy
1692 IF (PRESENT(calculate_forces)) c_forces = calculate_forces
1693
1694 IF (c_forces) THEN
1695 CALL timeset(routinen//'_forces', handle)
1696 ELSE
1697 CALL timeset(routinen, handle)
1698 END IF
1699
1700 cpassert(ASSOCIATED(qs_env))
1701
1702 CALL get_qs_env(qs_env, &
1703 ks_env=ks_env, &
1704 rho_changed=rho_changed, &
1705 s_mstruct_changed=s_mstruct_changed, &
1706 potential_changed=potential_changed, &
1707 forces_up_to_date=forces_up_to_date)
1708
1709 do_rebuild = .false.
1710 do_rebuild = do_rebuild .OR. rho_changed
1711 do_rebuild = do_rebuild .OR. s_mstruct_changed
1712 do_rebuild = do_rebuild .OR. potential_changed
1713 do_rebuild = do_rebuild .OR. (c_forces .AND. .NOT. forces_up_to_date)
1714
1715 IF (do_rebuild) THEN
1716 CALL evaluate_core_matrix_traces(qs_env)
1717
1718 ! the ks matrix will be rebuilt so this is fine now
1719 CALL set_ks_env(ks_env, potential_changed=.false.)
1720
1721 CALL rebuild_ks_matrix(qs_env, &
1722 calculate_forces=c_forces, &
1723 just_energy=energy_only, &
1724 print_active=print_active)
1725
1726 IF (.NOT. energy_only) THEN
1727 CALL set_ks_env(ks_env, &
1728 rho_changed=.false., &
1729 s_mstruct_changed=.false., &
1730 forces_up_to_date=forces_up_to_date .OR. c_forces)
1731 END IF
1732 END IF
1733
1734 CALL timestop(handle)
1735
1736 END SUBROUTINE qs_ks_update_qs_env
1737
1738! **************************************************************************************************
1739!> \brief Calculates the traces of the core matrices and the density matrix.
1740!> \param qs_env ...
1741!> \param rho_ao_ext ...
1742!> \author Ole Schuett
1743! **************************************************************************************************
1744 SUBROUTINE evaluate_core_matrix_traces(qs_env, rho_ao_ext)
1745 TYPE(qs_environment_type), POINTER :: qs_env
1746 TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
1747 POINTER :: rho_ao_ext
1748
1749 CHARACTER(LEN=*), PARAMETER :: routinen = 'evaluate_core_matrix_traces'
1750
1751 INTEGER :: handle
1752 REAL(kind=dp) :: energy_core_im
1753 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_t, rho_ao_kp
1754 TYPE(dft_control_type), POINTER :: dft_control
1755 TYPE(qs_energy_type), POINTER :: energy
1756 TYPE(qs_rho_type), POINTER :: rho
1757
1758 CALL timeset(routinen, handle)
1759 NULLIFY (energy, rho, dft_control, rho_ao_kp, matrixkp_t, matrixkp_h)
1760
1761 CALL get_qs_env(qs_env, &
1762 rho=rho, &
1763 energy=energy, &
1764 dft_control=dft_control, &
1765 kinetic_kp=matrixkp_t, &
1766 matrix_h_kp=matrixkp_h)
1767
1768 IF (PRESENT(rho_ao_ext)) THEN
1769 rho_ao_kp => rho_ao_ext
1770 ELSE
1771 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1772 END IF
1773
1774 CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy%core, dft_control%nspins)
1775
1776 ! Add the imaginary part in the RTP case
1777 IF (qs_env%run_rtp) THEN
1778 IF (dft_control%rtp_control%velocity_gauge) THEN
1779 CALL get_qs_env(qs_env, matrix_h_im_kp=matrixkp_h)
1780 CALL qs_rho_get(rho, rho_ao_im_kp=rho_ao_kp)
1781 CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy_core_im, dft_control%nspins)
1782 energy%core = energy%core - energy_core_im
1783 END IF
1784 END IF
1785
1786 ! kinetic energy
1787 IF (ASSOCIATED(matrixkp_t)) THEN
1788 CALL calculate_ptrace(matrixkp_t, rho_ao_kp, energy%kinetic, dft_control%nspins)
1789 END IF
1790
1791 CALL timestop(handle)
1792 END SUBROUTINE evaluate_core_matrix_traces
1793
1794! **************************************************************************************************
1795!> \brief Constructs a new Khon-Sham matrix
1796!> \param qs_env ...
1797!> \param calculate_forces ...
1798!> \param just_energy ...
1799!> \param print_active ...
1800!> \author Ole Schuett
1801! **************************************************************************************************
1802 SUBROUTINE rebuild_ks_matrix(qs_env, calculate_forces, just_energy, print_active)
1803 TYPE(qs_environment_type), POINTER :: qs_env
1804 LOGICAL, INTENT(IN) :: calculate_forces, just_energy
1805 LOGICAL, INTENT(IN), OPTIONAL :: print_active
1806
1807 CHARACTER(LEN=*), PARAMETER :: routinen = 'rebuild_ks_matrix'
1808
1809 INTEGER :: handle
1810 TYPE(dft_control_type), POINTER :: dft_control
1811 TYPE(qs_scf_env_type), POINTER :: scf_env
1812
1813 CALL timeset(routinen, handle)
1814 NULLIFY (dft_control, scf_env)
1815
1816 CALL get_qs_env(qs_env, dft_control=dft_control, scf_env=scf_env)
1817 IF (.NOT. just_energy .AND. ASSOCIATED(scf_env)) scf_env%ot_ks_valid = .false.
1818
1819 IF (dft_control%qs_control%semi_empirical) THEN
1820 CALL build_se_fock_matrix(qs_env, &
1821 calculate_forces=calculate_forces, &
1822 just_energy=just_energy)
1823
1824 ELSE IF (dft_control%qs_control%dftb) THEN
1825 CALL build_dftb_ks_matrix(qs_env, &
1826 calculate_forces=calculate_forces, &
1827 just_energy=just_energy)
1828
1829 ELSE IF (dft_control%qs_control%xtb) THEN
1830 IF (dft_control%qs_control%xtb_control%do_tblite) THEN
1831 CALL build_tblite_ks_matrix(qs_env, &
1832 calculate_forces=calculate_forces, &
1833 just_energy=just_energy)
1834 ELSE
1835 CALL build_xtb_ks_matrix(qs_env, &
1836 calculate_forces=calculate_forces, &
1837 just_energy=just_energy)
1838 END IF
1839 ELSE
1840 CALL qs_ks_build_kohn_sham_matrix(qs_env, &
1841 calculate_forces=calculate_forces, &
1842 just_energy=just_energy, &
1843 print_active=print_active)
1844 END IF
1845
1846 CALL timestop(handle)
1847
1848 END SUBROUTINE rebuild_ks_matrix
1849
1850! **************************************************************************************************
1851!> \brief Allocate ks_matrix if necessary, take current overlap matrix as template
1852!> \param qs_env ...
1853!> \param is_complex ...
1854!> \par History
1855!> refactoring 04.03.2011 [MI]
1856!> \author
1857! **************************************************************************************************
1858
1859 SUBROUTINE qs_ks_allocate_basics(qs_env, is_complex)
1860 TYPE(qs_environment_type), POINTER :: qs_env
1861 LOGICAL, INTENT(in) :: is_complex
1862
1863 CHARACTER(LEN=default_string_length) :: headline
1864 INTEGER :: ic, ispin, nimages, nspins
1865 LOGICAL :: do_kpoints
1866 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrixkp_im_ks, matrixkp_ks
1867 TYPE(dbcsr_type), POINTER :: refmatrix
1868 TYPE(dft_control_type), POINTER :: dft_control
1869 TYPE(kpoint_type), POINTER :: kpoints
1870 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1871 POINTER :: sab_orb
1872 TYPE(qs_ks_env_type), POINTER :: ks_env
1873
1874 NULLIFY (dft_control, ks_env, matrix_s_kp, sab_orb, matrixkp_ks, refmatrix, matrixkp_im_ks, kpoints)
1875
1876 CALL get_qs_env(qs_env, &
1877 dft_control=dft_control, &
1878 matrix_s_kp=matrix_s_kp, &
1879 ks_env=ks_env, &
1880 kpoints=kpoints, &
1881 do_kpoints=do_kpoints, &
1882 matrix_ks_kp=matrixkp_ks, &
1883 matrix_ks_im_kp=matrixkp_im_ks)
1884
1885 IF (do_kpoints) THEN
1886 CALL get_kpoint_info(kpoints, sab_nl=sab_orb)
1887 ELSE
1888 CALL get_qs_env(qs_env, sab_orb=sab_orb)
1889 END IF
1890
1891 nspins = dft_control%nspins
1892 nimages = dft_control%nimages
1893
1894 IF (.NOT. ASSOCIATED(matrixkp_ks)) THEN
1895 CALL dbcsr_allocate_matrix_set(matrixkp_ks, nspins, nimages)
1896 refmatrix => matrix_s_kp(1, 1)%matrix
1897 DO ispin = 1, nspins
1898 DO ic = 1, nimages
1899 IF (nspins > 1) THEN
1900 IF (ispin == 1) THEN
1901 headline = "KOHN-SHAM MATRIX FOR ALPHA SPIN"
1902 ELSE
1903 headline = "KOHN-SHAM MATRIX FOR BETA SPIN"
1904 END IF
1905 ELSE
1906 headline = "KOHN-SHAM MATRIX"
1907 END IF
1908 ALLOCATE (matrixkp_ks(ispin, ic)%matrix)
1909 CALL dbcsr_create(matrix=matrixkp_ks(ispin, ic)%matrix, template=refmatrix, &
1910 name=trim(headline), matrix_type=dbcsr_type_symmetric)
1911 CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_ks(ispin, ic)%matrix, sab_orb)
1912 CALL dbcsr_set(matrixkp_ks(ispin, ic)%matrix, 0.0_dp)
1913 END DO
1914 END DO
1915 CALL set_ks_env(ks_env, matrix_ks_kp=matrixkp_ks)
1916 END IF
1917
1918 IF (is_complex) THEN
1919 IF (.NOT. ASSOCIATED(matrixkp_im_ks)) THEN
1920 cpassert(nspins == SIZE(matrixkp_ks, 1))
1921 cpassert(nimages == SIZE(matrixkp_ks, 2))
1922 CALL dbcsr_allocate_matrix_set(matrixkp_im_ks, nspins, nimages)
1923 DO ispin = 1, nspins
1924 DO ic = 1, nimages
1925 IF (nspins > 1) THEN
1926 IF (ispin == 1) THEN
1927 headline = "IMAGINARY KOHN-SHAM MATRIX FOR ALPHA SPIN"
1928 ELSE
1929 headline = "IMAGINARY KOHN-SHAM MATRIX FOR BETA SPIN"
1930 END IF
1931 ELSE
1932 headline = "IMAGINARY KOHN-SHAM MATRIX"
1933 END IF
1934 ALLOCATE (matrixkp_im_ks(ispin, ic)%matrix)
1935 refmatrix => matrixkp_ks(ispin, ic)%matrix ! base on real part, but anti-symmetric
1936 CALL dbcsr_create(matrix=matrixkp_im_ks(ispin, ic)%matrix, template=refmatrix, &
1937 name=trim(headline), matrix_type=dbcsr_type_antisymmetric)
1938 CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_im_ks(ispin, ic)%matrix, sab_orb)
1939 CALL dbcsr_set(matrixkp_im_ks(ispin, ic)%matrix, 0.0_dp)
1940 END DO
1941 END DO
1942 CALL set_ks_env(ks_env, matrix_ks_im_kp=matrixkp_im_ks)
1943 END IF
1944 END IF
1945
1946 END SUBROUTINE qs_ks_allocate_basics
1947
1948END MODULE qs_ks_methods
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:96
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integral_ab
Definition pw_methods.F:222
pw_methods::pw_integrate_function
Definition pw_methods.F:116
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
pw_poisson_methods::pw_poisson_solve
Definition pw_poisson_methods.F:78
qs_core_energies::calculate_ptrace
Definition qs_core_energies.F:63
accint_weights_forces
Definition accint_weights_forces.F:12
accint_weights_forces::accint_weight_force
subroutine, public accint_weight_force(qs_env, rho, rho1, order, xc_section, triplet, force_scale)
...
Definition accint_weights_forces.F:93
admm_dm_methods
Contains ADMM methods which only require the density matrix.
Definition admm_dm_methods.F:14
admm_dm_methods::admm_dm_merge_ks_matrix
subroutine, public admm_dm_merge_ks_matrix(qs_env)
Entry methods: Merges auxiliary Kohn-Sham matrix into primary one.
Definition admm_dm_methods.F:94
admm_dm_methods::admm_dm_calc_rho_aux
subroutine, public admm_dm_calc_rho_aux(qs_env)
Entry methods: Calculates auxiliary density matrix from primary one.
Definition admm_dm_methods.F:57
admm_methods
Contains ADMM methods which require molecular orbitals.
Definition admm_methods.F:15
admm_methods::admm_mo_calc_rho_aux_kp
subroutine, public admm_mo_calc_rho_aux_kp(qs_env)
...
Definition admm_methods.F:284
admm_methods::admm_mo_merge_ks_matrix
subroutine, public admm_mo_merge_ks_matrix(qs_env)
...
Definition admm_methods.F:729
admm_methods::admm_update_ks_atom
subroutine, public admm_update_ks_atom(qs_env, calculate_forces)
Adds the GAPW exchange contribution to the aux_fit ks matrices.
Definition admm_methods.F:666
admm_methods::calc_admm_ovlp_forces_kp
subroutine, public calc_admm_ovlp_forces_kp(qs_env)
Calculate the forces due to the AUX/ORB basis overlap in ADMM, in the KP case.
Definition admm_methods.F:2731
admm_methods::admm_mo_calc_rho_aux
subroutine, public admm_mo_calc_rho_aux(qs_env)
...
Definition admm_methods.F:148
admm_methods::calc_admm_ovlp_forces
subroutine, public calc_admm_ovlp_forces(qs_env)
Calculate the forces due to the AUX/ORB basis overlap in ADMM.
Definition admm_methods.F:2690
admm_methods::calc_admm_mo_derivatives
subroutine, public calc_admm_mo_derivatives(qs_env, mo_derivs)
Calculate the derivative of the AUX_FIT mo, based on the ORB mo_derivs.
Definition admm_methods.F:2645
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:599
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:216
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_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_multiply
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
Definition cp_dbcsr_api.F:1096
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_filter
subroutine, public dbcsr_filter(matrix, eps)
...
Definition cp_dbcsr_api.F:658
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::dbcsr_copy_columns_hack
subroutine, public dbcsr_copy_columns_hack(matrix_b, matrix_a, ncol, source_start, target_start, para_env, blacs_env)
hack for dbcsr_copy_columns
Definition cp_dbcsr_operations.F:763
cp_ddapc
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition cp_ddapc.F:15
cp_ddapc::qs_ks_ddapc
subroutine, public qs_ks_ddapc(qs_env, auxbas_pw_pool, rho_tot_gspace, v_hartree_gspace, v_spin_ddapc_rest_r, energy, calculate_forces, ks_matrix, just_energy)
Set of methods using DDAPC charges.
Definition cp_ddapc.F:81
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
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_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
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_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:104
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:330
dft_plus_u
Add the DFT+U contribution to the Hamiltonian matrix.
Definition dft_plus_u.F:18
dft_plus_u::plus_u
subroutine, public plus_u(qs_env, matrix_h, matrix_w)
Add the DFT+U contribution to the Hamiltonian matrix. Wrapper routine for all "+U" methods.
Definition dft_plus_u.F:134
gapw_gspace_reconstruction
GAPW reciprocal-space reconstruction and its discrete adjoint.
Definition gapw_gspace_reconstruction.F:12
gapw_gspace_reconstruction::integrate_gapw_composite_vxc_gspace
subroutine, public integrate_gapw_composite_vxc_gspace(qs_env, auxbas_pw_pool, vxc_rho, vxc_tau, atom_force, strain_virial, one_center_contraction, one_center_rho_contraction, one_center_tau_contraction)
Apply the discrete adjoint of the common-grid GAPW hard-minus-soft reconstruction.
Definition gapw_gspace_reconstruction.F:382
gce_methods
Definition gce_methods.F:12
gce_methods::planar_averaged_v_hartree_3d
subroutine, public planar_averaged_v_hartree_3d(v_rspace, dft_control, do_gce, ref_esp, para_env)
calculate the planar averaged real space potential (e.g. Hartree potential) along the surface normal ...
Definition gce_methods.F:82
gce_methods::planar_counter_charge
subroutine, public planar_counter_charge(rho_tot_gspace, pcc_env, auxbas_pw_pool)
add the planar counter charge density to the total charge density
Definition gce_methods.F:190
hartree_local_methods
Definition hartree_local_methods.F:8
hartree_local_methods::vh_1c_gg_integrals
subroutine, public vh_1c_gg_integrals(qs_env, energy_hartree_1c, ecoul_1c, local_rho_set, para_env, tddft, local_rho_set_2nd, core_2nd)
Calculates one center GAPW Hartree energies and matrix elements Hartree potentials are input Takes po...
Definition hartree_local_methods.F:212
hartree_local_types
Definition hartree_local_types.F:8
hfx_ace_methods
Adaptively Compressed Exchange (ACE) operator for HFX. Reference: Lin, J. Chem. Theory Comput....
Definition hfx_ace_methods.F:43
hfx_ace_methods::hfx_ace_ks_matrix
subroutine, public hfx_ace_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, just_energy, v_rspace_new, v_tau_rspace, ace_rebuild_frequency, ext_xc_section)
Main ACE entry point, replacing hfx_ks_matrix in qs_ks_methods.
Definition hfx_ace_methods.F:151
hfx_admm_utils
Utilities for hfx and admm methods.
Definition hfx_admm_utils.F:17
hfx_admm_utils::hfx_admm_init
subroutine, public hfx_admm_init(qs_env, calculate_forces, ext_xc_section)
...
Definition hfx_admm_utils.F:165
hfx_admm_utils::hfx_ks_matrix
subroutine, public hfx_ks_matrix(qs_env, matrix_ks, rho, energy, calculate_forces, just_energy, v_rspace_new, v_tau_rspace, ext_xc_section)
Add the hfx contributions to the Hamiltonian.
Definition hfx_admm_utils.F:1143
hfx_admm_utils::hfx_ks_matrix_kp
subroutine, public hfx_ks_matrix_kp(qs_env, matrix_ks, energy, calculate_forces)
Add the HFX K-point contribution to the real-space Hamiltonians.
Definition hfx_admm_utils.F:965
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::smeagol_runtype_emtransport
integer, parameter, public smeagol_runtype_emtransport
Definition input_constants.F:1484
input_constants::do_ppl_grid
integer, parameter, public do_ppl_grid
Definition input_constants.F:294
input_constants::outer_scf_becke_constraint
integer, parameter, public outer_scf_becke_constraint
Definition input_constants.F:773
input_constants::outer_scf_hirshfeld_constraint
integer, parameter, public outer_scf_hirshfeld_constraint
Definition input_constants.F:773
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:735
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:708
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:1052
kg_correction
Routines for a Kim-Gordon-like partitioning into molecular subunits.
Definition kg_correction.F:14
kg_correction::kg_ekin_subset
subroutine, public kg_ekin_subset(qs_env, ks_matrix, ekin_mol, calc_force, do_kernel, pmat_ext)
Calculates the subsystem Hohenberg-Kohn kinetic energy and the forces.
Definition kg_correction.F:91
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, inversion_symmetry_only, symmetry_backend, symmetry_reduction_method, gamma_centered, lattice_fft)
Retrieve information from a kpoint environment.
Definition kpoint_types.F:613
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_energy
subroutine, public v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl_ri)
...
Definition lri_environment_methods.F:1221
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
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::abnormal_value
logical function, public abnormal_value(a)
determines if a value is not normal (e.g. for Inf and Nan) based on IO to work also under optimizatio...
Definition mathlib.F:159
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
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:67
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
qmmm_image_charge
Routines for image charge calculation within QM/MM.
Definition qmmm_image_charge.F:14
qmmm_image_charge::calculate_image_pot
subroutine, public calculate_image_pot(v_hartree_rspace, rho_hartree_gspace, energy, qmmm_env, qs_env)
determines coefficients by solving image_matrix*coeff=-pot_const by Gaussian elimination or in an ite...
Definition qmmm_image_charge.F:97
qmmm_image_charge::integrate_potential_devga_rspace
subroutine, public integrate_potential_devga_rspace(potential, coeff, forces, qmmm_env, qs_env)
calculates the image forces on the MM atoms
Definition qmmm_image_charge.F:454
qmmm_image_charge::add_image_pot_to_hartree_pot
subroutine, public add_image_pot_to_hartree_pot(v_hartree, v_metal, qs_env)
Add potential of metal (image charge pot) to Hartree Potential.
Definition qmmm_image_charge.F:921
qs_cdft_types
Defines CDFT control structures.
Definition qs_cdft_types.F:14
qs_cdft_utils
Utility subroutines for CDFT calculations.
Definition qs_cdft_utils.F:14
qs_cdft_utils::cdft_constraint_print
subroutine, public cdft_constraint_print(qs_env, electronic_charge)
Prints information about CDFT constraints.
Definition qs_cdft_utils.F:1117
qs_charges_types
container for information about total charges on the grids
Definition qs_charges_types.F:14
qs_core_energies
Calculation of the energies concerning the core charge distribution.
Definition qs_core_energies.F:14
qs_dftb_matrices
Calculation of Overlap and Hamiltonian matrices in DFTB.
Definition qs_dftb_matrices.F:12
qs_dftb_matrices::build_dftb_ks_matrix
subroutine, public build_dftb_ks_matrix(qs_env, calculate_forces, just_energy)
...
Definition qs_dftb_matrices.F:496
qs_efield_berry
Calculates the energy contribution and the mo_derivative of a static periodic electric field.
Definition qs_efield_berry.F:15
qs_efield_berry::qs_efield_berry_phase
subroutine, public qs_efield_berry_phase(qs_env, just_energy, calculate_forces)
...
Definition qs_efield_berry.F:107
qs_efield_local
Calculates the energy contribution and the mo_derivative of a static electric field (nonperiodic).
Definition qs_efield_local.F:15
qs_efield_local::qs_efield_local_operator
subroutine, public qs_efield_local_operator(qs_env, just_energy, calculate_forces)
...
Definition qs_efield_local.F:79
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, mimic, 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_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, 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, matrix_vhxc, 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, xcint_weights, 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, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_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, harris_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, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_force_types
Definition qs_force_types.F:13
qs_gapw_densities
Definition qs_gapw_densities.F:9
qs_gapw_densities::prepare_gapw_den
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external, pw_env_sub)
...
Definition qs_gapw_densities.F:57
qs_harris_types
Types needed for a for a Harris model calculation.
Definition qs_harris_types.F:14
qs_harris_utils
Harris method environment setup and handling.
Definition qs_harris_utils.F:14
qs_harris_utils::harris_set_potentials
subroutine, public harris_set_potentials(harris_env, vh_rspace, vxc_rspace)
...
Definition qs_harris_utils.F:580
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_ks_apply_restraints
Set of routines to apply restraints to the KS hamiltonian.
Definition qs_ks_apply_restraints.F:11
qs_ks_apply_restraints::qs_ks_s2_restraint
subroutine, public qs_ks_s2_restraint(dft_control, qs_env, matrix_s, energy, calculate_forces, just_energy)
...
Definition qs_ks_apply_restraints.F:187
qs_ks_apply_restraints::qs_ks_mulliken_restraint
subroutine, public qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, ks_matrix, matrix_s, rho, mulliken_order_p)
...
Definition qs_ks_apply_restraints.F:141
qs_ks_apply_restraints::qs_ks_cdft_constraint
subroutine, public qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
Apply a CDFT constraint.
Definition qs_ks_apply_restraints.F:59
qs_ks_atom
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
Definition qs_ks_atom.F:12
qs_ks_atom::update_ks_atom
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
Definition qs_ks_atom.F:110
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_methods.F:22
qs_ks_methods::rebuild_ks_matrix
subroutine, public rebuild_ks_matrix(qs_env, calculate_forces, just_energy, print_active)
Constructs a new Khon-Sham matrix.
Definition qs_ks_methods.F:1803
qs_ks_methods::evaluate_core_matrix_traces
subroutine, public evaluate_core_matrix_traces(qs_env, rho_ao_ext)
Calculates the traces of the core matrices and the density matrix.
Definition qs_ks_methods.F:1745
qs_ks_methods::qs_ks_update_qs_env
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
Definition qs_ks_methods.F:1674
qs_ks_methods::qs_ks_allocate_basics
subroutine, public qs_ks_allocate_basics(qs_env, is_complex)
Allocate ks_matrix if necessary, take current overlap matrix as template.
Definition qs_ks_methods.F:1860
qs_ks_methods::calc_rho_tot_gspace
subroutine, public calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
...
Definition qs_ks_methods.F:1530
qs_ks_methods::qs_ks_build_kohn_sham_matrix
subroutine, public qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, print_active, ext_ks_matrix, ext_xc_section)
routine where the real calculations are made: the KS matrix is calculated
Definition qs_ks_methods.F:206
qs_ks_qmmm_methods
Definition qs_ks_qmmm_methods.F:11
qs_ks_qmmm_methods::qmmm_calculate_energy
subroutine, public qmmm_calculate_energy(qs_env, rho, v_qmmm, qmmm_energy)
Computes the contribution to the total energy of the QM/MM electrostatic coupling.
Definition qs_ks_qmmm_methods.F:145
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:205
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, 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_vhxc, 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, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, 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_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:591
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:1109
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:2164
qs_ks_utils::compute_matrix_vxc_kp
subroutine, public compute_matrix_vxc_kp(qs_env, v_rspace, matrix_vxc_kp, gapw_full_basis)
Build the XC potential matrix for k-point/image-resolved KS matrices.
Definition qs_ks_utils.F:1490
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:168
qs_ks_utils::check_sum_energies
subroutine, public check_sum_energies(energy)
Check each term of the energy and sum to total.
Definition qs_ks_utils.F:1042
qs_ks_utils::compute_matrix_vxc
subroutine, public compute_matrix_vxc(qs_env, v_rspace, matrix_vxc, gapw_full_basis)
compute matrix_vxc, defined via the potential created by qs_vxc_create ignores things like tau functi...
Definition qs_ks_utils.F:1433
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:1573
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:1257
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:2051
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:923
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:627
qs_local_rho_types
Definition qs_local_rho_types.F:8
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, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:440
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_rho0_ggrid
Definition qs_rho0_ggrid.F:8
qs_rho0_ggrid::integrate_vhg0_rspace
subroutine, public integrate_vhg0_rspace(qs_env, v_rspace, para_env, calculate_forces, local_rho_set, local_rho_set_2nd, atener, kforce, my_pools, my_rs_descs)
...
Definition qs_rho0_ggrid.F:344
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
qs_sccs
Self-consistent continuum solvation (SCCS) model implementation.
Definition qs_sccs.F:29
qs_sccs::sccs
subroutine, public sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs, h_stress)
Self-consistent continuum solvation (SCCS) model implementation.
Definition qs_sccs.F:134
qs_sccs::sccs_core_fill_coefficients
subroutine, public sccs_core_fill_coefficients(qs_env, calpha, ccore, output_unit)
Exponents and coefficients of the SCCS core fill Gaussians per atomic kind, such that the fill densit...
Definition qs_sccs.F:1569
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
qs_vxc_atom
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition qs_vxc_atom.F:12
qs_vxc_atom::gapw_cdft_one_center
subroutine, public gapw_cdft_one_center(qs_env, energy_only, calculate_forces, values, electronic_charge, operator_group, rho_atom_operator_set)
Add the GAPW one-center correction to CDFT values and operators.
Definition qs_vxc_atom.F:2638
qs_vxc_atom::calculate_vxc_atom
subroutine, public calculate_vxc_atom(qs_env, energy_only, exc1, adiabatic_rescale_factor, kind_set_external, rho_atom_set_external, xc_section_external, calculate_forces, composite_vxc_rho, composite_vxc_tau, composite_reference_active, direct_valence_atom_grid, atom_composite_grid)
...
Definition qs_vxc_atom.F:231
qs_vxc
Definition qs_vxc.F:16
qs_vxc::qs_vxc_create
subroutine, public qs_vxc_create(ks_env, rho_struct, xc_section, vxc_rho, vxc_tau, exc, just_energy, edisp, dispersion_env, adiabatic_rescale_factor, pw_env_external, native_skala_atom_force, qs_env_external, native_gapw_composite_override, native_skala_defer_to_atom_composite)
calculates and allocates the xc potential, already reducing it to the dependence on rho and the one o...
Definition qs_vxc.F:120
rtp_admm_methods
Utilities for rtp in combination with admm methods adapted routines from admm_method (author Manuel G...
Definition rtp_admm_methods.F:15
rtp_admm_methods::rtp_admm_merge_ks_matrix
subroutine, public rtp_admm_merge_ks_matrix(qs_env)
...
Definition rtp_admm_methods.F:473
rtp_admm_methods::rtp_admm_calc_rho_aux
subroutine, public rtp_admm_calc_rho_aux(qs_env)
Compute the ADMM density matrix in case of rtp (complex MO's).
Definition rtp_admm_methods.F:85
se_fock_matrix
Calculation of the Fock matrix for SE methods.
Definition se_fock_matrix.F:17
se_fock_matrix::build_se_fock_matrix
subroutine, public build_se_fock_matrix(qs_env, calculate_forces, just_energy)
Construction of the Fock matrix for NDDO methods.
Definition se_fock_matrix.F:81
skala_gpw_functional
Experimental CP2K-native GPW real-space-grid path for SKALA TorchScript models.
Definition skala_gpw_functional.F:11
skala_gpw_functional::ensure_native_skala_grid_scope
subroutine, public ensure_native_skala_grid_scope(xc_section)
Enforce the currently implemented native SKALA GPW input scope.
Definition skala_gpw_functional.F:361
skala_gpw_functional::native_skala_gapw_composite_direct_ao
logical function, public native_skala_gapw_composite_direct_ao(xc_section)
Return true if the GAPW composite reference uses direct full-ORB collocation.
Definition skala_gpw_functional.F:260
skala_gpw_functional::xc_section_uses_gauxc_model
logical function, public xc_section_uses_gauxc_model(xc_section)
Return true if the GAUXC subsection requests a model evaluation.
Definition skala_gpw_functional.F:281
skala_gpw_functional::native_skala_uses_atom_composite_grid
logical function, public native_skala_uses_atom_composite_grid(xc_section)
Return true when native Skala uses atom-centered grids.
Definition skala_gpw_functional.F:230
skala_gpw_functional::skala_gapw_representation
integer function, public skala_gapw_representation(xc_section)
Return the pseudopotential GAPW representation selected for an active model.
Definition skala_gpw_functional.F:137
skala_gpw_functional::get_gauxc_section
type(section_vals_type) function, pointer, public get_gauxc_section(xc_section)
Return the first GAUXC functional subsection, if present.
Definition skala_gpw_functional.F:3061
skala_gpw_functional::native_skala_gapw_atom_composite_requested
logical function, public native_skala_gapw_atom_composite_requested(xc_section)
Return true when the explicit atom-centered composite reference is requested.
Definition skala_gpw_functional.F:205
skala_gpw_functional::native_skala_gapw_composite_reference
logical function, public native_skala_gapw_composite_reference(xc_section)
Return true if native SKALA should use the full GAPW ORB density on one common grid.
Definition skala_gpw_functional.F:181
skala_gpw_functional::xc_section_uses_native_skala_grid
logical function, public xc_section_uses_native_skala_grid(xc_section)
Return true if the GAUXC subsection requests the CP2K-native GPW grid path.
Definition skala_gpw_functional.F:118
smeagol_interface
CP2K+SMEAGOL interface.
Definition smeagol_interface.F:14
smeagol_interface::smeagol_shift_v_hartree
subroutine, public smeagol_shift_v_hartree(v_hartree_rspace, cell, hartreeleadsleft, hartreeleadsright, hartreeleadsbottom, vbias, zleft, zright, isexplicit_zright, isexplicit_bottom)
Align Hatree potential of semi-infinite leads to match bulk-transport calculation and apply external ...
Definition smeagol_interface.F:526
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
surface_dipole
Definition surface_dipole.F:9
surface_dipole::calc_dipsurf_potential
subroutine, public calc_dipsurf_potential(qs_env, energy)
compute the surface dipole and the correction to the hartree potential
Definition surface_dipole.F:59
tblite_ks_matrix
tblite matrix build
Definition tblite_ks_matrix.F:14
tblite_ks_matrix::build_tblite_ks_matrix
subroutine, public build_tblite_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
...
Definition tblite_ks_matrix.F:56
virial_types
Definition virial_types.F:13
xc_gauxc_functional
Definition xc_gauxc_functional.F:12
xc_gauxc_functional::gauxc_gapw_has_paw_pseudopotentials
logical function, public gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
Return whether GauXC GAPW mode sees pseudopotential one-center GAPW kinds.
Definition xc_gauxc_functional.F:439
xc_gauxc_functional::gauxc_gapw_all_all_electron
logical function, public gauxc_gapw_all_all_electron(qs_kind_set)
Return whether every GAPW kind uses an all-electron potential.
Definition xc_gauxc_functional.F:103
xc_gauxc_functional::gauxc_gapw_all_pseudopotentials
logical function, public gauxc_gapw_all_pseudopotentials(qs_kind_set)
Return whether every GAPW kind uses a pseudopotential.
Definition xc_gauxc_functional.F:410
xc_gauxc_functional::apply_gauxc
subroutine, public apply_gauxc(qs_env, xc_section, calculate_forces, rho_ao_external, exc_external, matrix_vxc_external)
Evaluate molecular GauXC, optionally for an explicit configuration density.
Definition xc_gauxc_functional.F:1022
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::skala_gapw_direct_valence
integer, parameter, public skala_gapw_direct_valence
Definition xc_input_constants.F:29
xc_input_constants::skala_gapw_paw_one_center
integer, parameter, public skala_gapw_paw_one_center
Definition xc_input_constants.F:29
xtb_ks_matrix
Calculation of KS matrix in xTB Reference: Stefan Grimme, Christoph Bannwarth, Philip Shushkov JCTC 1...
Definition xtb_ks_matrix.F:15
xtb_ks_matrix::build_xtb_ks_matrix
subroutine, public build_xtb_ks_matrix(qs_env, calculate_forces, just_energy, ext_ks_matrix)
...
Definition xtb_ks_matrix.F:73
admm_types::admm_type
stores some data used in wavefunction fitting
Definition admm_types.F:120
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:786
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_dbcsr_api::dbcsr_type
Definition cp_dbcsr_api.F:176
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
hartree_local_types::ecoul_1center_type
Definition hartree_local_types.F:25
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
kpoint_types::kpoint_type
Contains information about kpoints.
Definition kpoint_types.F:196
lri_environment_types::lri_density_type
Definition lri_environment_types.F:359
lri_environment_types::lri_environment_type
Definition lri_environment_types.F:271
lri_environment_types::lri_kind_type
Definition lri_environment_types.F:332
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_poisson_types::pw_poisson_type
environment for the poisson solver
Definition pw_poisson_types.F:123
pw_pool_types::pw_pool_type
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:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_cdft_types::cdft_control_type
Definition qs_cdft_types.F:218
qs_charges_types::qs_charges_type
Container for information about total charges on the grids.
Definition qs_charges_types.F:41
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:223
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_harris_types::harris_type
Contains information on the Harris method.
Definition qs_harris_types.F:66
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_local_rho_types::local_rho_type
Definition qs_local_rho_types.F:45
qs_mo_types::mo_set_type
Definition qs_mo_types.F:53
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:103
virial_types::virial_type
Definition virial_types.F:26