Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief routines that build the Kohn-Sham matrix (i.e calculate the coulomb
10 : !> and xc parts
11 : !> \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 : ! **************************************************************************************************
22 : MODULE 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 :
177 : CONTAINS
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 138581 : 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 138581 : 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 138581 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cdft_gapw_values
234 138581 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cdft_charge_correction, &
235 138581 : cdft_electronic_charge, &
236 138581 : 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 138581 : 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 138581 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_vhxc, matrix_vxc, mo_derivs
245 138581 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, ks_matrix_im, matrix_h, &
246 138581 : matrix_h_im, matrix_s, matrix_vxc_kp, &
247 138581 : my_rho, rho_ao
248 : TYPE(dft_control_type), POINTER :: dft_control
249 138581 : 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 138581 : 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 138581 : POINTER :: sab_orb
259 138581 : 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 415743 : TYPE(pw_r3d_rs_type), DIMENSION(2) :: v_hxc_rspace
266 277162 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace_embed, v_rspace_new, &
267 138581 : v_rspace_new_aux_fit, v_tau_rspace, &
268 138581 : 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 138581 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
273 138581 : 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 138581 : CALL timeset(routineN, handle)
283 138581 : NULLIFY (admm_env, atomic_kind_set, cell, dft_control, force, logger, mo_derivs, my_rho, &
284 138581 : rho_struct, para_env, pw_env, virial, vppl_rspace, &
285 138581 : ace_section, &
286 138581 : adiabatic_rescaling_section, hfx_sections, input, scf_section, &
287 138581 : xc_section, gauxc_section, matrix_h, matrix_h_im, matrix_s, auxbas_pw_pool, poisson_env, &
288 138581 : v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, matrix_vxc, &
289 138581 : matrix_vxc_kp, matrix_vhxc, &
290 138581 : vee, rho_nlcc, ks_env, ks_matrix, ks_matrix_im, rho, energy, rho_xc, rho_r, rho_ao, &
291 138581 : rho_core, sab_orb, particle_set, qs_kind_set, kpoints)
292 :
293 138581 : CPASSERT(ASSOCIATED(qs_env))
294 :
295 138581 : logger => cp_get_default_logger()
296 138581 : my_print = .TRUE.
297 138581 : IF (PRESENT(print_active)) my_print = print_active
298 138581 : use_gauxc_matrix = .FALSE.
299 138581 : native_gapw_atom_composite_reference = .FALSE.
300 138581 : native_gapw_atom_composite_requested = .FALSE.
301 138581 : native_gapw_all_electron_atom_composite = .FALSE.
302 138581 : native_direct_valence_atom_composite = .FALSE.
303 138581 : gapw_all_all_electron = .FALSE.
304 138581 : gapw_all_pseudopotentials = .FALSE.
305 138581 : gapw_mixed_core_representations = .FALSE.
306 138581 : native_gapw_composite_direct_ao = .FALSE.
307 138581 : native_gapw_composite_reference = .FALSE.
308 138581 : native_gapw_periodic_atom_composite = .FALSE.
309 138581 : native_grid_diagnostics = .FALSE.
310 138581 : native_skala_composite_contraction = 0.0_dp
311 138581 : native_skala_composite_rho_contraction = 0.0_dp
312 138581 : native_skala_composite_tau_contraction = 0.0_dp
313 138581 : 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 138581 : natom=natom)
346 :
347 : ! Reject force/stress requests before assembling even the baseline KS matrix.
348 138581 : IF (dft_control%low_spin_roks .AND. dft_control%qs_control%gapw) THEN
349 136 : IF (.NOT. virial%pv_numer .AND. &
350 : (virial%pv_calculate .OR. (virial%pv_availability .AND. calculate_forces))) THEN
351 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS stress is not implemented.")
352 : END IF
353 136 : IF (calculate_forces) THEN
354 0 : CALL cp_abort(__LOCATION__, "GAPW low spin ROKS forces are not implemented.")
355 : END IF
356 : END IF
357 :
358 138581 : CALL qs_rho_get(rho, rho_r=rho_r, rho_ao_kp=rho_ao)
359 :
360 138581 : nimages = dft_control%nimages
361 138581 : nspins = dft_control%nspins
362 :
363 : ! remap pointer to allow for non-kpoint external ks matrix
364 138581 : IF (PRESENT(ext_ks_matrix)) ks_matrix(1:nspins, 1:1) => ext_ks_matrix(1:nspins)
365 :
366 138581 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
367 :
368 138581 : adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
369 138581 : CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
370 138581 : just_energy_xc = just_energy
371 138581 : 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 44 : just_energy_xc = .TRUE.
375 : END IF
376 :
377 138581 : CPASSERT(ASSOCIATED(matrix_h))
378 138581 : CPASSERT(ASSOCIATED(matrix_s))
379 138581 : CPASSERT(ASSOCIATED(rho))
380 138581 : CPASSERT(ASSOCIATED(pw_env))
381 138581 : CPASSERT(SIZE(ks_matrix, 1) > 0)
382 138581 : 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 138581 : qs_env%cp_ddapc_ewald%do_solvation
389 :
390 : ! Check if LRIGPW is used
391 138581 : lrigpw = dft_control%qs_control%lrigpw
392 138581 : rigpw = dft_control%qs_control%rigpw
393 138581 : IF (rigpw) THEN
394 26 : CPASSERT(nimages == 1)
395 : END IF
396 26 : IF (lrigpw .AND. rigpw) THEN
397 0 : CPABORT(" LRI and RI are not compatible")
398 : END IF
399 :
400 : ! Check for GAPW method : additional terms for local densities
401 138581 : gapw = dft_control%qs_control%gapw
402 138581 : gapw_xc = dft_control%qs_control%gapw_xc
403 138581 : IF (gapw_xc .AND. gapw) THEN
404 0 : CPABORT(" GAPW and GAPW_XC are not compatible")
405 : END IF
406 138581 : IF ((gapw .AND. lrigpw) .OR. (gapw_xc .AND. lrigpw)) THEN
407 0 : CPABORT(" GAPW/GAPW_XC and LRIGPW are not compatible")
408 : END IF
409 138581 : IF ((gapw .AND. rigpw) .OR. (gapw_xc .AND. rigpw)) THEN
410 0 : CPABORT(" GAPW/GAPW_XC and RIGPW are not compatible")
411 : END IF
412 :
413 138581 : do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
414 138581 : IF (do_ppl) THEN
415 60 : CPASSERT(.NOT. gapw)
416 60 : CALL get_qs_env(qs_env=qs_env, vppl=vppl_rspace)
417 : END IF
418 :
419 138581 : IF (gapw_xc) THEN
420 4446 : CPASSERT(ASSOCIATED(rho_xc))
421 : END IF
422 :
423 : ! gets the tmp grids
424 138581 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
425 :
426 138581 : IF (gapw .AND. (poisson_env%parameters%solver == pw_poisson_implicit)) THEN
427 0 : CPABORT("The implicit Poisson solver cannot be used in conjunction with GAPW.")
428 : END IF
429 :
430 : ! *** Prepare densities for gapw ***
431 138581 : IF (gapw .OR. gapw_xc) THEN
432 28908 : CALL prepare_gapw_den(qs_env, do_rho0=(.NOT. gapw_xc))
433 : END IF
434 :
435 : ! Calculate the Hartree potential
436 138581 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
437 138581 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
438 :
439 138581 : 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 138581 : (.NOT. gapw) .AND. (.NOT. gapw_xc) .AND. &
444 : (.NOT. (poisson_env%parameters%solver == pw_poisson_implicit))) THEN
445 936 : CALL pw_zero(rho_tot_gspace)
446 936 : 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 936 : v_hartree_gspace)
449 936 : CALL pw_zero(rho_tot_gspace)
450 936 : CALL pw_zero(v_hartree_gspace)
451 : END IF
452 :
453 : ! Get the total density in g-space [ions + electrons]
454 138581 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
455 :
456 138581 : IF (qs_env%scf_control%gce%do_gce .AND. .NOT. dft_control%do_pcc) THEN
457 0 : CPABORT("GCE requires DFT%PLANAR_COUNTER_CHARGE to define the countercharge plane.")
458 : END IF
459 :
460 : ! Add the planar counter charge density
461 138581 : IF (dft_control%do_pcc) THEN
462 412 : CALL planar_counter_charge(rho_tot_gspace, dft_control%pcc_control, auxbas_pw_pool)
463 : END IF
464 :
465 138581 : IF (my_print) THEN
466 138525 : CALL print_densities(qs_env, rho)
467 : END IF
468 :
469 138581 : IF (dft_control%do_sccs) THEN
470 : ! Self-consistent continuum solvation (SCCS) model
471 : NULLIFY (v_sccs_rspace)
472 262 : ALLOCATE (v_sccs_rspace)
473 262 : CALL auxbas_pw_pool%create_pw(v_sccs_rspace)
474 :
475 262 : IF (poisson_env%parameters%solver == pw_poisson_implicit) THEN
476 0 : CPABORT("The implicit Poisson solver cannot be used together with SCCS.")
477 : END IF
478 :
479 262 : IF (use_virial .AND. calculate_forces) THEN
480 : CALL sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs_rspace, &
481 0 : h_stress=h_stress)
482 0 : virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
483 0 : virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
484 : ELSE
485 262 : 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 138319 : IF (use_virial .AND. calculate_forces) THEN
491 544 : 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 544 : rho_core=rho_core)
495 7072 : virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
496 7072 : 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 137775 : v_hartree_gspace, rho_core=rho_core)
500 : END IF
501 : END IF
502 :
503 138581 : IF (dft_control%do_paep .OR. qs_env%scf_control%gce%do_gce) THEN
504 388 : 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 388 : 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 138581 : 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 2340 : just_energy)
514 : ELSE
515 136241 : dft_control%qs_control%ddapc_explicit_potential = .FALSE.
516 136241 : dft_control%qs_control%ddapc_restraint_is_spin = .FALSE.
517 136241 : IF (.NOT. just_energy) THEN
518 125747 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
519 125747 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
520 : END IF
521 : END IF
522 138581 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
523 :
524 138581 : IF (dft_control%correct_surf_dip) THEN
525 110 : IF (dft_control%surf_dip_correct_switch) THEN
526 110 : CALL calc_dipsurf_potential(qs_env, energy)
527 110 : 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 138581 : just_energy, calculate_forces, auxbas_pw_pool)
534 :
535 : ! Check if CDFT constraint is needed
536 138581 : CALL qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
537 :
538 : ! Adds the External Potential if requested
539 138581 : IF (dft_control%apply_external_potential) THEN
540 : ! Compute the energy due to the external potential
541 : ee_ener = 0.0_dp
542 728 : DO ispin = 1, nspins
543 728 : ee_ener = ee_ener + pw_integral_ab(rho_r(ispin), vee)
544 : END DO
545 364 : IF (.NOT. just_energy) THEN
546 364 : IF (gapw) THEN
547 : CALL get_qs_env(qs_env=qs_env, &
548 : rho0_s_rs=rho0_s_rs, &
549 42 : rhoz_cneo_s_rs=rhoz_cneo_s_rs)
550 42 : CPASSERT(ASSOCIATED(rho0_s_rs))
551 42 : IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
552 0 : CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs)
553 : END IF
554 42 : ee_ener = ee_ener + pw_integral_ab(rho0_s_rs, vee)
555 42 : IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
556 0 : 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 364 : energy%ee = -ee_ener
562 : END IF
563 :
564 : ! Adds the QM/MM potential
565 138581 : 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 6334 : qmmm_energy=energy%qmmm_el)
570 6334 : 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 60 : qs_env=qs_env)
576 60 : 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 60 : qs_env=qs_env)
580 60 : 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 20 : qmmm_env=qs_env%qmmm_env_qm, qs_env=qs_env)
585 : END IF
586 : END IF
587 60 : CALL qs_env%ks_qmmm_env%v_metal_rspace%release()
588 60 : DEALLOCATE (qs_env%ks_qmmm_env%v_metal_rspace)
589 : END IF
590 6334 : IF (.NOT. just_energy) THEN
591 : CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, &
592 6244 : v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, scale=1.0_dp)
593 : END IF
594 : END IF
595 138581 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
596 :
597 : ! SMEAGOL interface
598 138581 : IF (dft_control%smeagol_control%smeagol_enabled .AND. &
599 : dft_control%smeagol_control%run_type == smeagol_runtype_emtransport) THEN
600 0 : 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 0 : dft_control%smeagol_control%aux%isexplicit_HartreeLeadsBottom)
610 : END IF
611 :
612 : ! calculate the density matrix for the fitted mo_coeffs
613 138581 : IF (dft_control%do_admm) THEN
614 13710 : IF (PRESENT(ext_xc_section)) THEN
615 0 : CALL hfx_admm_init(qs_env, calculate_forces, ext_xc_section)
616 : ELSE
617 13710 : CALL hfx_admm_init(qs_env, calculate_forces)
618 : END IF
619 :
620 13710 : IF (dft_control%do_admm_mo) THEN
621 13496 : IF (qs_env%run_rtp) THEN
622 114 : CALL rtp_admm_calc_rho_aux(qs_env)
623 : ELSE
624 13382 : IF (dokp) THEN
625 242 : CALL admm_mo_calc_rho_aux_kp(qs_env)
626 : ELSE
627 13140 : CALL admm_mo_calc_rho_aux(qs_env)
628 : END IF
629 : END IF
630 214 : ELSE IF (dft_control%do_admm_dm) THEN
631 214 : CALL admm_dm_calc_rho_aux(qs_env)
632 : END IF
633 : END IF
634 :
635 : ! only activate stress calculation if
636 138581 : 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 138581 : IF (dft_control%do_admm) THEN
642 13710 : CALL get_qs_env(qs_env, admm_env=admm_env)
643 13710 : xc_section => admm_env%xc_section_aux
644 13710 : 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 13710 : just_energy=just_energy_xc)
650 :
651 13710 : 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 4970 : calculate_forces=calculate_forces)
658 :
659 : END IF
660 :
661 13710 : NULLIFY (rho_struct)
662 :
663 13710 : IF (use_virial .AND. calculate_forces) THEN
664 22 : vscale = 1.0_dp
665 : !Note: ADMMS and ADMMP stress tensor only for closed-shell calculations
666 22 : IF (admm_env%do_admms) vscale = admm_env%gsi(1)**(2.0_dp/3.0_dp)
667 22 : IF (admm_env%do_admmp) vscale = admm_env%gsi(1)**2
668 286 : virial%pv_exc = virial%pv_exc - vscale*virial%pv_xc
669 286 : 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 13710 : xc_section => admm_env%xc_section_primary
673 : ELSE
674 124871 : 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 124871 : IF (PRESENT(ext_xc_section)) xc_section => ext_xc_section
677 : END IF
678 :
679 138581 : IF (gapw_xc) THEN
680 4446 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
681 : ELSE
682 134135 : CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
683 : END IF
684 :
685 : ! zmp
686 138581 : IF (dft_control%apply_external_density .OR. dft_control%apply_external_vxc) THEN
687 106 : energy%exc = 0.0_dp
688 106 : CALL calculate_zmp_potential(qs_env, v_rspace_new, rho, exc=energy%exc)
689 : ELSE
690 : ! Embedding potential (runs regardless of XC method)
691 138475 : IF (dft_control%apply_embed_pot) THEN
692 822 : NULLIFY (v_rspace_embed)
693 822 : energy%embed_corr = 0.0_dp
694 : CALL get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, &
695 822 : energy%embed_corr, just_energy)
696 : END IF
697 :
698 : ! Everything else, either via GauXC or manual XC computation
699 138475 : native_skala_evaluator = xc_section_uses_native_skala_grid(xc_section)
700 138475 : paw_composite_representation = .FALSE.
701 138475 : 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 310 : gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
705 : native_skala_evaluator = native_skala_evaluator .OR. &
706 310 : paw_composite_representation
707 : END IF
708 138475 : IF (dft_control%use_gauxc) THEN
709 906 : IF (native_skala_evaluator) THEN
710 300 : IF (gapw .AND. dft_control%low_spin_roks) THEN
711 : CALL cp_abort(__LOCATION__, &
712 0 : "GAPW low spin ROKS requires the molecular GauXC grid, not native Skala grids.")
713 : END IF
714 300 : CALL ensure_native_skala_grid_scope(xc_section)
715 300 : gapw_representation = skala_gapw_representation(xc_section)
716 300 : IF (gapw .OR. gapw_xc) THEN
717 168 : gapw_all_all_electron = gauxc_gapw_all_all_electron(qs_kind_set)
718 168 : gapw_all_pseudopotentials = gauxc_gapw_all_pseudopotentials(qs_kind_set)
719 : gapw_mixed_core_representations = &
720 168 : .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 300 : gapw_all_pseudopotentials))
728 : native_gapw_composite_reference = &
729 300 : 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 300 : (gapw .OR. gapw_xc)
733 : native_gapw_atom_composite_reference = &
734 300 : native_gapw_atom_composite_requested .OR. paw_composite_representation
735 : native_gapw_periodic_atom_composite = &
736 300 : native_skala_uses_atom_composite_grid(xc_section)
737 300 : 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 0 : "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 300 : 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 300 : native_gapw_periodic_atom_composite
753 : native_gapw_atom_composite_reference = &
754 : native_gapw_atom_composite_reference .OR. &
755 300 : native_gapw_all_electron_atom_composite
756 300 : IF (paw_composite_representation) THEN
757 136 : IF (ANY(cell%perd /= 0)) THEN
758 : native_gapw_composite_reference = &
759 104 : .NOT. native_gapw_periodic_atom_composite
760 : native_gapw_atom_composite_reference = &
761 104 : native_gapw_periodic_atom_composite
762 : ELSE
763 8 : native_gapw_composite_reference = .FALSE.
764 8 : native_gapw_atom_composite_reference = .TRUE.
765 : END IF
766 : END IF
767 300 : IF (native_gapw_composite_reference .AND. native_gapw_atom_composite_reference) THEN
768 : CALL cp_abort(__LOCATION__, &
769 0 : "Select only one native-grid GAPW composite reference route.")
770 : END IF
771 : native_gapw_composite_direct_ao = native_gapw_composite_reference .AND. &
772 300 : native_skala_gapw_composite_direct_ao(xc_section)
773 300 : gauxc_section => get_gauxc_section(xc_section)
774 300 : IF (ASSOCIATED(gauxc_section)) THEN
775 : CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_DIAGNOSTICS", &
776 300 : l_val=native_grid_diagnostics)
777 : END IF
778 300 : IF ((.NOT. do_kpoints) .AND. nimages /= 1) THEN
779 : CALL cp_abort(__LOCATION__, &
780 : "Native SKALA grid evaluation supports multiple images only "// &
781 0 : "for k-point calculations.")
782 : END IF
783 300 : IF (do_kpoints) THEN
784 36 : CPASSERT(ASSOCIATED(kpoints))
785 36 : gauxc_section => get_gauxc_section(xc_section)
786 36 : CPASSERT(ASSOCIATED(gauxc_section))
787 36 : CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_USE_CUDA", l_val=native_grid_use_cuda)
788 36 : IF (.NOT. native_grid_use_cuda) THEN
789 36 : 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 3 : "runtime stack, or select NATIVE_GRID_USE_CUDA T.")
794 3 : native_grid_cpu_kpoints_warned = .TRUE.
795 : END IF
796 : END IF
797 : END IF
798 300 : IF (dft_control%roks) THEN
799 0 : CPABORT("Native SKALA grid evaluation does not support ROKS.")
800 : END IF
801 300 : IF (dft_control%do_admm) THEN
802 0 : 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 300 : native_skala_restore_exc = calculate_forces
808 300 : IF (native_skala_restore_exc) THEN
809 66 : native_skala_exc_scf = energy%exc
810 66 : native_skala_total_scf = energy%total
811 : END IF
812 300 : IF (calculate_forces) THEN
813 198 : ALLOCATE (native_skala_atom_force(3, natom))
814 66 : native_skala_atom_force = 0.0_dp
815 66 : 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 58 : 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 8 : 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 66 : ASSOCIATED(v_rspace_new) .AND. ASSOCIATED(v_tau_rspace)) THEN
834 2 : IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
835 0 : output_unit = cp_logger_get_default_io_unit()
836 0 : IF (output_unit > 0) THEN
837 0 : DO iatom = 1, natom
838 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,I0,3(1X,ES20.12))") &
839 0 : "SKALA_GPW| Native regular-grid atom force", iatom, &
840 0 : native_skala_atom_force(:, iatom)
841 : END DO
842 : END IF
843 : END IF
844 2 : IF (use_virial) THEN
845 2 : 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 2 : one_center_tau_contraction=native_skala_composite_tau_contraction)
853 2 : IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
854 0 : output_unit = cp_logger_get_default_io_unit()
855 0 : IF (output_unit > 0) THEN
856 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
857 0 : "SKALA_GPW| Composite-direction one-center virial"
858 0 : DO iatom = 1, 3
859 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,3ES20.10)") &
860 0 : "SKALA_GPW|", native_skala_composite_virial(iatom, :)
861 : END DO
862 : END IF
863 : END IF
864 26 : 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 0 : one_center_tau_contraction=native_skala_composite_tau_contraction)
872 : END IF
873 : END IF
874 66 : IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
875 7 : output_unit = cp_logger_get_default_io_unit()
876 7 : IF (output_unit > 0) THEN
877 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
878 7 : "SKALA_GPW| Composite-direction one-center VXC contraction", &
879 14 : native_skala_composite_contraction
880 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
881 7 : "SKALA_GPW| Composite-direction one-center rho contraction", &
882 14 : native_skala_composite_rho_contraction
883 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
884 7 : "SKALA_GPW| Composite-direction one-center tau contraction", &
885 14 : native_skala_composite_tau_contraction
886 21 : DO iatom = 1, natom
887 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,I0,3(1X,ES20.12))") &
888 21 : "SKALA_GPW| Native atom force", iatom, native_skala_atom_force(:, iatom)
889 : END DO
890 : END IF
891 : END IF
892 66 : CPASSERT(ASSOCIATED(force))
893 66 : CPASSERT(ASSOCIATED(atomic_kind_set))
894 66 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
895 198 : DO iatom = 1, natom
896 132 : ikind = kind_of(iatom)
897 132 : atom_a = atom_of_kind(iatom)
898 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + &
899 594 : native_skala_atom_force(:, iatom)
900 : END DO
901 66 : 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 234 : native_direct_valence_atom_composite))
911 : IF (native_gapw_composite_reference .AND. .NOT. native_gapw_composite_direct_ao .AND. &
912 234 : 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 2 : one_center_tau_contraction=native_skala_composite_tau_contraction)
918 2 : IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
919 0 : output_unit = cp_logger_get_default_io_unit()
920 0 : IF (output_unit > 0) THEN
921 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
922 0 : "SKALA_GPW| Composite-direction one-center VXC contraction", &
923 0 : native_skala_composite_contraction
924 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
925 0 : "SKALA_GPW| Composite-direction one-center rho contraction", &
926 0 : native_skala_composite_rho_contraction
927 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
928 0 : "SKALA_GPW| Composite-direction one-center tau contraction", &
929 0 : native_skala_composite_tau_contraction
930 : END IF
931 : END IF
932 : END IF
933 : END IF
934 300 : IF (native_skala_restore_exc) energy%exc = native_skala_exc_scf
935 300 : IF (native_direct_valence_atom_composite) THEN
936 124 : CPASSERT(ASSOCIATED(v_rspace_new))
937 124 : 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 124 : composite_vxc_tau=v_tau_rspace, direct_valence_atom_grid=.TRUE.)
942 124 : energy%exc = 0.0_dp
943 176 : ELSE IF ((gapw .OR. gapw_xc) .AND. .NOT. native_gapw_composite_reference) THEN
944 144 : IF (native_gapw_atom_composite_reference) THEN
945 128 : CPASSERT(ASSOCIATED(v_rspace_new))
946 128 : 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 128 : composite_vxc_tau=v_tau_rspace, atom_composite_grid=.TRUE.)
951 128 : 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 16 : calculate_forces=calculate_forces)
956 : END IF
957 : END IF
958 300 : IF (edisp /= 0.0_dp) energy%dispersion = edisp
959 300 : IF (qs_env%requires_matrix_vxc .AND. ASSOCIATED(v_rspace_new)) THEN
960 0 : 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 0 : gapw_full_basis=native_gapw_composite_direct_ao)
964 0 : 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 0 : gapw_full_basis=native_gapw_composite_direct_ao)
969 0 : CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
970 : END IF
971 : END IF
972 : ELSE
973 606 : use_gauxc_matrix = .TRUE.
974 606 : CALL apply_gauxc(qs_env, xc_section, calculate_forces)
975 606 : 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 2 : 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 137569 : just_energy=just_energy_xc)
986 137569 : IF (edisp /= 0.0_dp) energy%dispersion = edisp
987 137569 : IF (qs_env%requires_matrix_vxc .AND. ASSOCIATED(v_rspace_new)) THEN
988 2 : CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace_new, matrix_vxc=matrix_vxc)
989 2 : CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
990 : END IF
991 :
992 137569 : 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 28598 : calculate_forces=calculate_forces)
996 : END IF
997 : END IF
998 : END IF
999 :
1000 29146 : IF (gapw .AND. dft_control%qs_control%cdft) THEN
1001 1440 : ALLOCATE (cdft_gapw_values(SIZE(cdft_control%group)))
1002 480 : IF (cdft_control%atomic_charges) THEN
1003 : ALLOCATE (cdft_charge_correction(natom, nspins), &
1004 532 : cdft_electronic_charge(cdft_control%natoms, nspins))
1005 : CALL gapw_cdft_one_center(qs_env, just_energy, calculate_forces, &
1006 76 : cdft_gapw_values, cdft_charge_correction)
1007 228 : DO ispin = 1, nspins
1008 624 : DO iatom = 1, cdft_control%natoms
1009 : cdft_electronic_charge(iatom, ispin) = &
1010 548 : pw_integral_ab(cdft_control%charge(iatom), rho_r(ispin), local_only=.TRUE.)
1011 : END DO
1012 : END DO
1013 76 : CALL para_env%sum(cdft_electronic_charge)
1014 228 : DO ispin = 1, nspins
1015 624 : DO iatom = 1, cdft_control%natoms
1016 : cdft_electronic_charge(iatom, ispin) = &
1017 : cdft_electronic_charge(iatom, ispin) + &
1018 548 : cdft_charge_correction(cdft_control%atoms(iatom), ispin)
1019 : END DO
1020 : END DO
1021 : ELSE
1022 404 : CALL gapw_cdft_one_center(qs_env, just_energy, calculate_forces, cdft_gapw_values)
1023 : END IF
1024 1040 : cdft_control%value = cdft_control%value + cdft_gapw_values
1025 1040 : energy%cdft = energy%cdft + DOT_PRODUCT(cdft_control%strength, cdft_gapw_values)
1026 480 : IF (cdft_control%atomic_charges) THEN
1027 76 : CALL cdft_constraint_print(qs_env, cdft_electronic_charge)
1028 76 : DEALLOCATE (cdft_charge_correction, cdft_electronic_charge)
1029 : ELSE
1030 404 : CALL cdft_constraint_print(qs_env, cdft_electronic_charge)
1031 : END IF
1032 480 : DEALLOCATE (cdft_gapw_values)
1033 : END IF
1034 :
1035 138581 : IF (dft_control%mtlr_dft_with_perturbation .AND. .NOT. just_energy) THEN
1036 1844 : IF (do_adiabatic_rescaling .OR. use_gauxc_matrix) THEN
1037 0 : CPABORT("MTLR requires an XC potential on the real-space grid.")
1038 : END IF
1039 1844 : IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
1040 80 : ALLOCATE (matrix_vhxc(nspins))
1041 48 : DO ispin = 1, nspins
1042 32 : NULLIFY (matrix_vhxc(ispin)%matrix)
1043 32 : ALLOCATE (matrix_vhxc(ispin)%matrix)
1044 : CALL dbcsr_create(matrix_vhxc(ispin)%matrix, &
1045 32 : template=ks_matrix(1, 1)%matrix)
1046 32 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_vhxc(ispin)%matrix, sab_orb)
1047 48 : CALL dbcsr_set(matrix_vhxc(ispin)%matrix, 0.0_dp)
1048 : END DO
1049 16 : CALL set_ks_env(ks_env, matrix_vhxc=matrix_vhxc)
1050 : ELSE
1051 1828 : CPASSERT(SIZE(matrix_vhxc) == nspins)
1052 5484 : DO ispin = 1, nspins
1053 5484 : CPASSERT(ASSOCIATED(matrix_vhxc(ispin)%matrix))
1054 : END DO
1055 : END IF
1056 5532 : DO ispin = 1, nspins
1057 3688 : CALL auxbas_pw_pool%create_pw(v_hxc_rspace(ispin))
1058 3688 : CALL pw_zero(v_hxc_rspace(ispin))
1059 3688 : CALL pw_axpy(v_hartree_rspace, v_hxc_rspace(ispin))
1060 3688 : IF (ASSOCIATED(v_rspace_new)) THEN
1061 : CALL pw_axpy(v_rspace_new(ispin), &
1062 : v_hxc_rspace(ispin), &
1063 3688 : v_rspace_new(ispin)%pw_grid%dvol)
1064 : END IF
1065 3688 : 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 3688 : calculate_forces=.FALSE.)
1070 5532 : 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 138581 : IF (qs_env%harris_method) THEN
1076 116 : CALL get_qs_env(qs_env, harris_env=harris_env)
1077 116 : CALL harris_set_potentials(harris_env, v_hartree_rspace, v_rspace_new)
1078 : END IF
1079 :
1080 138581 : NULLIFY (rho_struct)
1081 138581 : IF (use_virial .AND. calculate_forces) THEN
1082 7072 : virial%pv_exc = virial%pv_exc - virial%pv_xc
1083 7072 : virial%pv_virial = virial%pv_virial - virial%pv_xc
1084 : END IF
1085 :
1086 : ! *** Add Hartree-Fock contribution if required ***
1087 138581 : hfx_sections => section_vals_get_subs_vals(xc_section, "HF")
1088 138581 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
1089 :
1090 138581 : ace_active = .FALSE.
1091 138581 : ace_rebuild_frequency = 1
1092 :
1093 138581 : IF (do_hfx) THEN
1094 31652 : ace_section => section_vals_get_subs_vals(hfx_sections, "ACE")
1095 31652 : IF (ASSOCIATED(ace_section)) THEN
1096 31652 : CALL section_vals_val_get(ace_section, "ACTIVE", l_val=ace_active)
1097 31652 : CALL section_vals_val_get(ace_section, "REBUILD_FREQUENCY", i_val=ace_rebuild_frequency)
1098 : END IF
1099 : END IF
1100 :
1101 138581 : IF (do_hfx) THEN
1102 31652 : IF (dokp) THEN
1103 364 : IF (ace_active) THEN
1104 0 : CPABORT("ACE-HFX for k-points is not implemented yet")
1105 : ELSE
1106 364 : 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 31288 : 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 48 : ace_rebuild_frequency, ext_xc_section=xc_section)
1115 : ELSE
1116 : CALL hfx_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, &
1117 31240 : 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 138581 : IF (do_ppl .AND. calculate_forces) THEN
1123 12 : CPASSERT(.NOT. gapw)
1124 26 : DO ispin = 1, nspins
1125 26 : 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 138581 : .NOT. native_gapw_atom_composite_reference .AND. &
1132 : .NOT. native_direct_valence_atom_composite) THEN
1133 92 : DO ispin = 1, nspins
1134 46 : CALL integrate_rho_nlcc(v_rspace_new(ispin), qs_env)
1135 92 : 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 138581 : IF (dft_control%qs_control%do_kg .AND. just_energy) THEN
1141 :
1142 12 : CPASSERT(nimages == 1)
1143 12 : ksmat => ks_matrix(:, 1)
1144 12 : CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces, do_kernel=.FALSE.)
1145 :
1146 : ! subtract kg corr from the total energy
1147 12 : energy%exc = energy%exc - ekin_mol
1148 :
1149 : END IF
1150 :
1151 : ! *** Single atom contributions ***
1152 138581 : IF (.NOT. just_energy) THEN
1153 127697 : IF (calculate_forces) THEN
1154 : ! Getting nuclear force contribution from the core charge density
1155 6029 : IF ((poisson_env%parameters%solver == pw_poisson_implicit) .AND. &
1156 : (poisson_env%parameters%dielectric_params%dielec_core_correction)) THEN
1157 28 : BLOCK
1158 : TYPE(pw_r3d_rs_type) :: v_minus_veps
1159 28 : CALL auxbas_pw_pool%create_pw(v_minus_veps)
1160 28 : CALL pw_copy(v_hartree_rspace, v_minus_veps)
1161 28 : CALL pw_axpy(poisson_env%implicit_env%v_eps, v_minus_veps, -v_hartree_rspace%pw_grid%dvol)
1162 28 : CALL integrate_v_core_rspace(v_minus_veps, qs_env)
1163 28 : CALL auxbas_pw_pool%give_back_pw(v_minus_veps)
1164 : END BLOCK
1165 : ELSE
1166 6001 : 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 6029 : IF (dft_control%do_sccs) THEN
1173 6 : IF (dft_control%sccs_control%core_fill) THEN
1174 : BLOCK
1175 2 : REAL(KIND=dp), DIMENSION(:), ALLOCATABLE :: calpha, ccore
1176 : TYPE(pw_r3d_rs_type) :: v_minus_vsccs
1177 2 : CALL sccs_core_fill_coefficients(qs_env, calpha, ccore)
1178 2 : CALL auxbas_pw_pool%create_pw(v_minus_vsccs)
1179 2 : CALL pw_copy(v_sccs_rspace, v_minus_vsccs)
1180 2 : CALL pw_scale(v_minus_vsccs, -1.0_dp)
1181 2 : CALL integrate_v_core_rspace(v_minus_vsccs, qs_env, calpha=calpha, ccore=ccore)
1182 2 : CALL auxbas_pw_pool%give_back_pw(v_minus_vsccs)
1183 4 : DEALLOCATE (calpha, ccore)
1184 : END BLOCK
1185 : END IF
1186 : END IF
1187 : END IF
1188 :
1189 127697 : 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 216975 : DO ispin = 1, nspins
1193 1421783 : DO img = 1, nimages
1194 1204808 : CALL dbcsr_get_info(ks_matrix(ispin, img)%matrix, name=name) ! keep the name
1195 1323298 : 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 98485 : IF (qs_env%run_rtp) THEN
1200 2088 : IF (dft_control%rtp_control%velocity_gauge) THEN
1201 182 : CPASSERT(ASSOCIATED(matrix_h_im))
1202 182 : CPASSERT(ASSOCIATED(ks_matrix_im))
1203 364 : DO ispin = 1, nspins
1204 546 : DO img = 1, nimages
1205 182 : CALL dbcsr_get_info(ks_matrix_im(ispin, img)%matrix, name=name) ! keep the name
1206 364 : 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 127697 : IF (use_virial .AND. calculate_forces) THEN
1214 7072 : 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 127697 : 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 127697 : cdft_control, calculate_forces)
1224 :
1225 127697 : IF (use_gauxc_matrix) THEN
1226 594 : IF (dokp) THEN
1227 0 : CALL get_qs_env(qs_env=qs_env, matrix_vxc_kp=matrix_vxc_kp)
1228 0 : CPASSERT(ASSOCIATED(matrix_vxc_kp))
1229 0 : DO ispin = 1, nspins
1230 0 : DO img = 1, nimages
1231 : CALL dbcsr_add(ks_matrix(ispin, img)%matrix, matrix_vxc_kp(ispin, img)%matrix, &
1232 0 : 1.0_dp, 1.0_dp)
1233 : END DO
1234 : END DO
1235 : ELSE
1236 594 : CALL get_qs_env(qs_env=qs_env, matrix_vxc=matrix_vxc)
1237 594 : CPASSERT(ASSOCIATED(matrix_vxc))
1238 594 : CPASSERT(nimages == 1)
1239 1230 : DO ispin = 1, nspins
1240 1230 : 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 127697 : IF (gapw .OR. gapw_xc) THEN
1246 27936 : IF (calculate_forces) THEN
1247 790 : IF (gapw_xc) THEN
1248 124 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
1249 : ELSE
1250 666 : CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
1251 : END IF
1252 790 : NULLIFY (rho1)
1253 790 : IF (dft_control%use_gauxc .AND. (gapw .OR. gapw_xc) .AND. &
1254 : .NOT. native_skala_evaluator) THEN
1255 2 : gauxc_model_none = .FALSE.
1256 2 : gauxc_section => get_gauxc_section(xc_section)
1257 2 : IF (ASSOCIATED(gauxc_section)) THEN
1258 2 : CALL section_vals_val_get(gauxc_section, "MODEL", c_val=gauxc_model_name)
1259 2 : gauxc_model_name = ADJUSTL(gauxc_model_name)
1260 2 : CALL uppercase(gauxc_model_name)
1261 : gauxc_model_none = (TRIM(gauxc_model_name) == "" .OR. &
1262 2 : TRIM(gauxc_model_name) == "NONE")
1263 : END IF
1264 2 : IF (gauxc_model_none .AND. &
1265 : (gapw_xc .OR. gauxc_gapw_has_paw_pseudopotentials(qs_kind_set))) THEN
1266 0 : CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section)
1267 : END IF
1268 : ELSE
1269 788 : CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section)
1270 : END IF
1271 : !
1272 790 : IF (dft_control%do_admm) THEN
1273 92 : CALL get_qs_env(qs_env, admm_env=admm_env)
1274 92 : xc_section => admm_env%xc_section_aux
1275 92 : CALL get_admm_env(admm_env, rho_aux_fit=rho_struct)
1276 92 : vscale = 1.0_dp
1277 92 : IF (admm_env%do_admmp) THEN
1278 8 : vscale = admm_env%gsi(1)**2
1279 84 : ELSE IF (admm_env%do_admms) THEN
1280 6 : vscale = admm_env%gsi(1)**(2.0_dp/3.0_dp)
1281 : END IF
1282 92 : 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 127697 : IF (use_virial .AND. calculate_forces) THEN
1288 7072 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1289 : END IF
1290 127697 : IF (dft_control%qs_control%do_kg) THEN
1291 900 : CPASSERT(nimages == 1)
1292 900 : ksmat => ks_matrix(:, 1)
1293 :
1294 900 : IF (use_virial .AND. calculate_forces) THEN
1295 208 : pv_loc = virial%pv_virial
1296 : END IF
1297 :
1298 900 : CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces, do_kernel=.FALSE.)
1299 : ! subtract kg corr from the total energy
1300 900 : energy%exc = energy%exc - ekin_mol
1301 :
1302 : ! virial corrections
1303 900 : IF (use_virial .AND. calculate_forces) THEN
1304 :
1305 : ! Integral contribution
1306 208 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1307 :
1308 : ! GGA contribution
1309 208 : virial%pv_exc = virial%pv_exc + virial%pv_xc
1310 208 : virial%pv_virial = virial%pv_virial + virial%pv_xc
1311 208 : 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 10884 : IF (ASSOCIATED(v_rspace_new)) THEN
1319 0 : DO ispin = 1, SIZE(v_rspace_new)
1320 0 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
1321 : END DO
1322 0 : DEALLOCATE (v_rspace_new)
1323 : END IF
1324 10884 : IF (ASSOCIATED(v_tau_rspace)) THEN
1325 0 : DO ispin = 1, SIZE(v_tau_rspace)
1326 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1327 : END DO
1328 0 : 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 138581 : IF (dft_control%qs_control%ddapc_explicit_potential) THEN
1337 164 : CALL auxbas_pw_pool%give_back_pw(v_spin_ddapc_rest_r)
1338 164 : DEALLOCATE (v_spin_ddapc_rest_r)
1339 : END IF
1340 :
1341 138581 : IF (calculate_forces .AND. dft_control%qs_control%cdft) THEN
1342 126 : IF (.NOT. cdft_control%transfer_pot) THEN
1343 230 : DO iatom = 1, SIZE(cdft_control%group)
1344 124 : CALL auxbas_pw_pool%give_back_pw(cdft_control%group(iatom)%weight)
1345 230 : DEALLOCATE (cdft_control%group(iatom)%weight)
1346 : END DO
1347 106 : IF (cdft_control%atomic_charges) THEN
1348 78 : DO iatom = 1, cdft_control%natoms
1349 78 : CALL auxbas_pw_pool%give_back_pw(cdft_control%charge(iatom))
1350 : END DO
1351 26 : DEALLOCATE (cdft_control%charge)
1352 : END IF
1353 106 : IF (cdft_control%type == outer_scf_becke_constraint .AND. &
1354 : cdft_control%becke_control%cavity_confine) THEN
1355 88 : IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN
1356 64 : CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
1357 : ELSE
1358 24 : DEALLOCATE (cdft_control%becke_control%cavity_mat)
1359 : END IF
1360 18 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1361 6 : IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) THEN
1362 0 : CALL auxbas_pw_pool%give_back_pw(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)
1363 : END IF
1364 : END IF
1365 106 : IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment)
1366 106 : cdft_control%save_pot = .FALSE.
1367 106 : cdft_control%need_pot = .TRUE.
1368 106 : cdft_control%external_control = .FALSE.
1369 : END IF
1370 : END IF
1371 :
1372 138581 : IF (dft_control%do_sccs) THEN
1373 262 : CALL auxbas_pw_pool%give_back_pw(v_sccs_rspace)
1374 262 : DEALLOCATE (v_sccs_rspace)
1375 : END IF
1376 :
1377 138581 : IF (gapw) THEN
1378 24462 : 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 42 : v_qmmm=vee, scale=-1.0_dp)
1382 : END IF
1383 24462 : 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 24462 : 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 24462 : core_2nd=.FALSE.)
1389 : ! CNEO quantum nuclear core energy (kinetic + Z*erfc(r)/r potential from classical nuclei)
1390 24462 : energy%core_cneo = 0.0_dp
1391 24462 : IF (ASSOCIATED(local_rho_set%rhoz_cneo_set)) THEN
1392 184 : DO iatom = 1, SIZE(local_rho_set%rhoz_cneo_set)
1393 184 : 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 138581 : IF (gapw .OR. gapw_xc) THEN
1399 : ! Single atom contributions in the KS matrix ***
1400 28908 : CALL update_ks_atom(qs_env, ks_matrix, rho_ao, calculate_forces)
1401 28908 : 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 4970 : 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 138581 : ks_matrix, matrix_s, rho, mulliken_order_p)
1411 :
1412 : ! Add DFT+U contribution, if requested
1413 138581 : IF (dft_control%dft_plus_u) THEN
1414 6324 : IF (just_energy) THEN
1415 746 : CALL plus_u(qs_env=qs_env)
1416 : ELSE
1417 5578 : CALL plus_u(qs_env=qs_env, matrix_h=ks_matrix)
1418 : END IF
1419 : ELSE
1420 132257 : 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 306859 : DO ispin = 1, nspins
1425 1803069 : DO img = 1, nimages
1426 : CALL dbcsr_filter(ks_matrix(ispin, img)%matrix, &
1427 1664488 : 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 138581 : IF (dft_control%do_admm_mo) THEN
1433 13496 : IF (qs_env%run_rtp) THEN
1434 114 : CALL rtp_admm_merge_ks_matrix(qs_env)
1435 : ELSE
1436 13382 : CALL admm_mo_merge_ks_matrix(qs_env)
1437 : END IF
1438 125085 : ELSE IF (dft_control%do_admm_dm) THEN
1439 214 : CALL admm_dm_merge_ks_matrix(qs_env)
1440 : END IF
1441 :
1442 : ! External field (nonperiodic case)
1443 138581 : 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 138581 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy .AND. .NOT. qs_env%run_rtp) THEN
1450 46925 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs)
1451 46925 : CPASSERT(nimages == 1)
1452 46925 : ksmat => ks_matrix(:, 1)
1453 46925 : CALL calc_mo_derivatives(qs_env, ksmat, mo_derivs)
1454 : END IF
1455 :
1456 : ! ADMM overlap forces
1457 138581 : IF (calculate_forces .AND. dft_control%do_admm) THEN
1458 320 : IF (dokp) THEN
1459 32 : CALL calc_admm_ovlp_forces_kp(qs_env)
1460 : ELSE
1461 288 : 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 138581 : 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 138581 : calculate_forces, auxbas_pw_pool)
1472 :
1473 : ! Periodic external field
1474 138581 : 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 138581 : energy, calculate_forces, just_energy)
1479 :
1480 138581 : IF (do_ppl) THEN
1481 : ! update core energy for grid based local pseudopotential
1482 60 : ecore_ppl = 0._dp
1483 126 : DO ispin = 1, nspins
1484 126 : ecore_ppl = ecore_ppl + pw_integral_ab(vppl_rspace, rho_r(ispin))
1485 : END DO
1486 60 : energy%core = energy%core + ecore_ppl
1487 : END IF
1488 :
1489 138581 : IF (lrigpw) THEN
1490 : ! update core energy for ppl_ri method
1491 486 : CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density)
1492 486 : IF (lri_env%ppl_ri) THEN
1493 8 : ecore_ppl = 0._dp
1494 16 : DO ispin = 1, nspins
1495 8 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1496 16 : CALL v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl)
1497 : END DO
1498 8 : energy%core = energy%core + ecore_ppl
1499 : END IF
1500 : END IF
1501 :
1502 : ! Sum all energy terms to obtain the total energy
1503 138581 : IF (native_skala_restore_exc) THEN
1504 66 : energy%total = native_skala_total_scf
1505 66 : IF (abnormal_value(energy%total)) THEN
1506 : CALL cp_abort(__LOCATION__, &
1507 0 : "SCF energy restored by native Skala is an abnormal value (NaN/Inf).")
1508 : END IF
1509 : ELSE
1510 138515 : CALL check_sum_energies(energy)
1511 : END IF
1512 :
1513 : ! Print detailed energy
1514 138581 : IF (my_print) THEN
1515 138525 : CALL print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
1516 : END IF
1517 :
1518 138581 : CALL timestop(handle)
1519 :
1520 415743 : 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 142387 : 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 142387 : 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 142387 : my_skip = .FALSE.
1543 950 : IF (PRESENT(skip_nuclear_density)) my_skip = skip_nuclear_density
1544 :
1545 142387 : CALL qs_rho_get(rho, rho_g=rho_g)
1546 142387 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
1547 :
1548 142387 : IF (.NOT. my_skip) THEN
1549 141447 : NULLIFY (rho_core)
1550 141447 : CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1551 141447 : IF (dft_control%qs_control%gapw) THEN
1552 24780 : NULLIFY (rho0_s_gs, rhoz_cneo_s_gs)
1553 24780 : CALL get_qs_env(qs_env=qs_env, rho0_s_gs=rho0_s_gs, rhoz_cneo_s_gs=rhoz_cneo_s_gs)
1554 24780 : CPASSERT(ASSOCIATED(rho0_s_gs))
1555 24780 : CALL pw_copy(rho0_s_gs, rho_tot_gspace)
1556 24780 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
1557 48 : CALL pw_axpy(rhoz_cneo_s_gs, rho_tot_gspace)
1558 : END IF
1559 24780 : IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
1560 2876 : CALL pw_axpy(rho_core, rho_tot_gspace)
1561 : END IF
1562 : ELSE
1563 116667 : CALL pw_copy(rho_core, rho_tot_gspace)
1564 : END IF
1565 312917 : DO ispin = 1, dft_control%nspins
1566 312917 : CALL pw_axpy(rho_g(ispin), rho_tot_gspace)
1567 : END DO
1568 141447 : CALL get_qs_env(qs_env=qs_env, qs_charges=qs_charges)
1569 141447 : qs_charges%total_rho_gspace = pw_integrate_function(rho_tot_gspace, isign=-1)
1570 : ELSE
1571 1884 : DO ispin = 1, dft_control%nspins
1572 1884 : CALL pw_axpy(rho_g(ispin), rho_tot_gspace)
1573 : END DO
1574 : END IF
1575 :
1576 142387 : 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 46925 : 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 46925 : 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 46925 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
1600 :
1601 46925 : 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 46925 : mos=mo_array)
1606 :
1607 104107 : DO ispin = 1, SIZE(mo_derivs)
1608 :
1609 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, &
1610 57182 : 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 57182 : 0.0_dp, mo_derivs(ispin)%matrix)
1613 :
1614 104107 : IF (dft_control%restricted) THEN
1615 : ! only the first mo_set are actual variables, but we still need both
1616 564 : CPASSERT(ispin == 1)
1617 564 : 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 564 : CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
1623 564 : CPASSERT(uniform_occupation)
1624 564 : CALL get_mo_set(mo_set=mo_array(2), uniform_occupation=uniform_occupation)
1625 564 : CPASSERT(uniform_occupation)
1626 :
1627 : ! The beta-spin might have fewer orbitals than alpa-spin...
1628 : ! create temporary matrices with beta_nmo columns
1629 564 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff_b)
1630 564 : CALL dbcsr_create(mo_derivs2_tmp1, template=mo_coeff_b)
1631 :
1632 : ! calculate beta derivatives
1633 564 : 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 564 : CALL dbcsr_create(mo_derivs2_tmp2, template=mo_derivs(1)%matrix)
1637 564 : 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 564 : blacs_env=mo_array(1)%mo_coeff%matrix_struct%context)
1644 :
1645 : ! add beta contribution to alpa mo_derivs
1646 564 : CALL dbcsr_add(mo_derivs(1)%matrix, mo_derivs2_tmp2, 1.0_dp, 1.0_dp)
1647 564 : CALL dbcsr_release(mo_derivs2_tmp1)
1648 564 : CALL dbcsr_release(mo_derivs2_tmp2)
1649 : END IF
1650 : END DO
1651 :
1652 46925 : IF (dft_control%do_admm_mo) THEN
1653 6864 : CALL calc_admm_mo_derivatives(qs_env, mo_derivs)
1654 : END IF
1655 :
1656 46925 : 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 284462 : 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 284462 : NULLIFY (ks_env)
1687 284462 : unit_nr = cp_logger_get_default_io_unit()
1688 :
1689 284462 : c_forces = .FALSE.
1690 284462 : energy_only = .FALSE.
1691 284462 : IF (PRESENT(just_energy)) energy_only = just_energy
1692 284462 : IF (PRESENT(calculate_forces)) c_forces = calculate_forces
1693 :
1694 284462 : IF (c_forces) THEN
1695 10619 : CALL timeset(routineN//'_forces', handle)
1696 : ELSE
1697 273843 : CALL timeset(routineN, handle)
1698 : END IF
1699 :
1700 284462 : 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 284462 : forces_up_to_date=forces_up_to_date)
1708 :
1709 284462 : do_rebuild = .FALSE.
1710 284462 : do_rebuild = do_rebuild .OR. rho_changed
1711 9123 : do_rebuild = do_rebuild .OR. s_mstruct_changed
1712 9123 : do_rebuild = do_rebuild .OR. potential_changed
1713 9123 : do_rebuild = do_rebuild .OR. (c_forces .AND. .NOT. forces_up_to_date)
1714 :
1715 : IF (do_rebuild) THEN
1716 275711 : CALL evaluate_core_matrix_traces(qs_env)
1717 :
1718 : ! the ks matrix will be rebuilt so this is fine now
1719 275711 : 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 275711 : print_active=print_active)
1725 :
1726 275711 : IF (.NOT. energy_only) THEN
1727 : CALL set_ks_env(ks_env, &
1728 : rho_changed=.FALSE., &
1729 : s_mstruct_changed=.FALSE., &
1730 508107 : forces_up_to_date=forces_up_to_date .OR. c_forces)
1731 : END IF
1732 : END IF
1733 :
1734 284462 : CALL timestop(handle)
1735 :
1736 284462 : 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 300607 : 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 300607 : 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 300607 : CALL timeset(routineN, handle)
1759 300607 : 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 300607 : matrix_h_kp=matrixkp_h)
1767 :
1768 300607 : IF (PRESENT(rho_ao_ext)) THEN
1769 24118 : rho_ao_kp => rho_ao_ext
1770 : ELSE
1771 276489 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1772 : END IF
1773 :
1774 300607 : CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy%core, dft_control%nspins)
1775 :
1776 : ! Add the imaginary part in the RTP case
1777 300607 : IF (qs_env%run_rtp) THEN
1778 3350 : IF (dft_control%rtp_control%velocity_gauge) THEN
1779 182 : CALL get_qs_env(qs_env, matrix_h_im_kp=matrixkp_h)
1780 182 : CALL qs_rho_get(rho, rho_ao_im_kp=rho_ao_kp)
1781 182 : CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy_core_im, dft_control%nspins)
1782 182 : energy%core = energy%core - energy_core_im
1783 : END IF
1784 : END IF
1785 :
1786 : ! kinetic energy
1787 300607 : IF (ASSOCIATED(matrixkp_t)) THEN
1788 138213 : CALL calculate_ptrace(matrixkp_t, rho_ao_kp, energy%kinetic, dft_control%nspins)
1789 : END IF
1790 :
1791 300607 : CALL timestop(handle)
1792 300607 : 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 275731 : 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 275731 : CALL timeset(routineN, handle)
1814 275731 : NULLIFY (dft_control, scf_env)
1815 :
1816 275731 : CALL get_qs_env(qs_env, dft_control=dft_control, scf_env=scf_env)
1817 275731 : IF (.NOT. just_energy .AND. ASSOCIATED(scf_env)) scf_env%ot_ks_valid = .FALSE.
1818 :
1819 275731 : IF (dft_control%qs_control%semi_empirical) THEN
1820 : CALL build_se_fock_matrix(qs_env, &
1821 : calculate_forces=calculate_forces, &
1822 41346 : just_energy=just_energy)
1823 :
1824 234385 : ELSE IF (dft_control%qs_control%dftb) THEN
1825 : CALL build_dftb_ks_matrix(qs_env, &
1826 : calculate_forces=calculate_forces, &
1827 31420 : just_energy=just_energy)
1828 :
1829 202965 : ELSE IF (dft_control%qs_control%xtb) THEN
1830 64754 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
1831 : CALL build_tblite_ks_matrix(qs_env, &
1832 : calculate_forces=calculate_forces, &
1833 27480 : just_energy=just_energy)
1834 : ELSE
1835 : CALL build_xtb_ks_matrix(qs_env, &
1836 : calculate_forces=calculate_forces, &
1837 37274 : 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 138211 : print_active=print_active)
1844 : END IF
1845 :
1846 275731 : CALL timestop(handle)
1847 :
1848 275731 : 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 29806 : 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 29806 : 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 29806 : POINTER :: sab_orb
1872 : TYPE(qs_ks_env_type), POINTER :: ks_env
1873 :
1874 29806 : 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 29806 : matrix_ks_im_kp=matrixkp_im_ks)
1884 :
1885 29806 : IF (do_kpoints) THEN
1886 4106 : CALL get_kpoint_info(kpoints, sab_nl=sab_orb)
1887 : ELSE
1888 25700 : CALL get_qs_env(qs_env, sab_orb=sab_orb)
1889 : END IF
1890 :
1891 29806 : nspins = dft_control%nspins
1892 29806 : nimages = dft_control%nimages
1893 :
1894 29806 : IF (.NOT. ASSOCIATED(matrixkp_ks)) THEN
1895 29756 : CALL dbcsr_allocate_matrix_set(matrixkp_ks, nspins, nimages)
1896 29756 : refmatrix => matrix_s_kp(1, 1)%matrix
1897 63786 : DO ispin = 1, nspins
1898 413878 : DO ic = 1, nimages
1899 350092 : IF (nspins > 1) THEN
1900 75444 : IF (ispin == 1) THEN
1901 37722 : headline = "KOHN-SHAM MATRIX FOR ALPHA SPIN"
1902 : ELSE
1903 37722 : headline = "KOHN-SHAM MATRIX FOR BETA SPIN"
1904 : END IF
1905 : ELSE
1906 274648 : headline = "KOHN-SHAM MATRIX"
1907 : END IF
1908 350092 : ALLOCATE (matrixkp_ks(ispin, ic)%matrix)
1909 : CALL dbcsr_create(matrix=matrixkp_ks(ispin, ic)%matrix, template=refmatrix, &
1910 350092 : name=TRIM(headline), matrix_type=dbcsr_type_symmetric)
1911 350092 : CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_ks(ispin, ic)%matrix, sab_orb)
1912 384122 : CALL dbcsr_set(matrixkp_ks(ispin, ic)%matrix, 0.0_dp)
1913 : END DO
1914 : END DO
1915 29756 : CALL set_ks_env(ks_env, matrix_ks_kp=matrixkp_ks)
1916 : END IF
1917 :
1918 29806 : IF (is_complex) THEN
1919 174 : IF (.NOT. ASSOCIATED(matrixkp_im_ks)) THEN
1920 174 : CPASSERT(nspins == SIZE(matrixkp_ks, 1))
1921 174 : CPASSERT(nimages == SIZE(matrixkp_ks, 2))
1922 174 : CALL dbcsr_allocate_matrix_set(matrixkp_im_ks, nspins, nimages)
1923 372 : DO ispin = 1, nspins
1924 570 : DO ic = 1, nimages
1925 198 : IF (nspins > 1) THEN
1926 48 : IF (ispin == 1) THEN
1927 24 : headline = "IMAGINARY KOHN-SHAM MATRIX FOR ALPHA SPIN"
1928 : ELSE
1929 24 : headline = "IMAGINARY KOHN-SHAM MATRIX FOR BETA SPIN"
1930 : END IF
1931 : ELSE
1932 150 : headline = "IMAGINARY KOHN-SHAM MATRIX"
1933 : END IF
1934 198 : ALLOCATE (matrixkp_im_ks(ispin, ic)%matrix)
1935 198 : 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 198 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric)
1938 198 : CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_im_ks(ispin, ic)%matrix, sab_orb)
1939 396 : CALL dbcsr_set(matrixkp_im_ks(ispin, ic)%matrix, 0.0_dp)
1940 : END DO
1941 : END DO
1942 174 : CALL set_ks_env(ks_env, matrix_ks_im_kp=matrixkp_im_ks)
1943 : END IF
1944 : END IF
1945 :
1946 29806 : END SUBROUTINE qs_ks_allocate_basics
1947 :
1948 : END MODULE qs_ks_methods
|
