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, compute_matrix_vxc, compute_matrix_vxc_kp, &
130 : get_embed_potential_energy, low_spin_roks, print_densities, print_detailed_energy, &
131 : 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 : USE qs_vxc, ONLY: qs_vxc_create
141 : USE qs_vxc_atom, ONLY: calculate_vxc_atom,&
142 : gapw_cdft_one_center
143 : USE rtp_admm_methods, ONLY: rtp_admm_calc_rho_aux,&
144 : rtp_admm_merge_ks_matrix
145 : USE se_fock_matrix, ONLY: build_se_fock_matrix
146 : USE skala_gpw_functional, ONLY: ensure_native_skala_grid_scope,&
147 : get_gauxc_section,&
148 : native_skala_gapw_atom_composite_requested,&
149 : native_skala_gapw_composite_direct_ao,&
150 : native_skala_gapw_composite_reference,&
151 : skala_gapw_representation,&
152 : 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_has_paw_pseudopotentials
160 : USE xc_input_constants, ONLY: skala_gapw_paw_one_center
161 : USE xtb_ks_matrix, ONLY: build_xtb_ks_matrix
162 : #include "./base/base_uses.f90"
163 :
164 : IMPLICIT NONE
165 :
166 : PRIVATE
167 :
168 : LOGICAL, PARAMETER :: debug_this_module = .TRUE.
169 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_methods'
170 :
171 : PUBLIC :: calc_rho_tot_gspace, qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, &
172 : qs_ks_allocate_basics, evaluate_core_matrix_traces, rebuild_ks_matrix
173 :
174 : CONTAINS
175 :
176 : ! **************************************************************************************************
177 : !> \brief routine where the real calculations are made: the
178 : !> KS matrix is calculated
179 : !> \param qs_env the qs_env to update
180 : !> \param calculate_forces if true calculate the quantities needed
181 : !> to calculate the forces. Defaults to false.
182 : !> \param just_energy if true updates the energies but not the
183 : !> ks matrix. Defaults to false
184 : !> \param print_active ...
185 : !> \param ext_ks_matrix ...
186 : !> \param ext_xc_section ...
187 : !> \par History
188 : !> 06.2002 moved from qs_scf to qs_ks_methods, use of ks_env
189 : !> new did_change scheme [fawzi]
190 : !> 10.2002 introduced pools, uses updated rho as input, LSD [fawzi]
191 : !> 10.2004 build_kohn_sham matrix now also computes the derivatives
192 : !> of the total energy wrt to the MO coefs, if instructed to
193 : !> do so. This appears useful for orbital dependent functionals
194 : !> where the KS matrix alone (however this might be defined)
195 : !> does not contain the info to construct this derivative.
196 : !> \author Matthias Krack
197 : !> \note
198 : !> make rho, energy and qs_charges optional, defaulting
199 : !> to qs_env components?
200 : ! **************************************************************************************************
201 124119 : SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
202 : print_active, ext_ks_matrix, ext_xc_section)
203 : TYPE(qs_environment_type), POINTER :: qs_env
204 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
205 : LOGICAL, INTENT(IN), OPTIONAL :: print_active
206 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
207 : POINTER :: ext_ks_matrix
208 : TYPE(section_vals_type), OPTIONAL, POINTER :: ext_xc_section
209 :
210 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_build_kohn_sham_matrix'
211 :
212 : CHARACTER(len=default_string_length) :: gauxc_model_name, name
213 : INTEGER :: ace_rebuild_frequency, atom_a, handle, &
214 : iatom, ikind, img, ispin, natom, &
215 : nimages, nspins, output_unit
216 124119 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
217 : LOGICAL :: ace_active, do_adiabatic_rescaling, do_ddapc, do_hfx, do_kpoints, do_ppl, dokp, &
218 : gapw, gapw_xc, gauxc_model_none, just_energy_xc, lrigpw, my_print, &
219 : native_gapw_atom_composite_reference, native_gapw_atom_composite_requested, &
220 : native_gapw_composite_direct_ao, native_gapw_composite_reference, &
221 : native_grid_diagnostics, native_grid_use_cuda, native_skala_evaluator, &
222 : native_skala_restore_exc, paw_composite_representation, rigpw, use_gauxc_matrix, &
223 : use_virial
224 : LOGICAL, SAVE :: native_grid_cpu_kpoints_warned = .FALSE.
225 : REAL(KIND=dp) :: ecore_ppl, edisp, ee_ener, ekin_mol, mulliken_order_p, &
226 : native_skala_composite_contraction, native_skala_composite_rho_contraction, &
227 : native_skala_composite_tau_contraction, native_skala_exc_scf, native_skala_total_scf, &
228 : vscale
229 124119 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cdft_gapw_values
230 124119 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cdft_charge_correction, &
231 124119 : cdft_electronic_charge, &
232 124119 : native_skala_atom_force
233 : REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, native_skala_composite_virial, &
234 : pv_loc
235 : TYPE(admm_type), POINTER :: admm_env
236 124119 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
237 : TYPE(cdft_control_type), POINTER :: cdft_control
238 : TYPE(cell_type), POINTER :: cell
239 : TYPE(cp_logger_type), POINTER :: logger
240 124119 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_vhxc, matrix_vxc, mo_derivs
241 124119 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, ks_matrix_im, matrix_h, &
242 124119 : matrix_h_im, matrix_s, matrix_vxc_kp, &
243 124119 : my_rho, rho_ao
244 : TYPE(dft_control_type), POINTER :: dft_control
245 124119 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
246 : TYPE(harris_type), POINTER :: harris_env
247 : TYPE(kpoint_type), POINTER :: kpoints
248 : TYPE(local_rho_type), POINTER :: local_rho_set
249 : TYPE(lri_density_type), POINTER :: lri_density
250 : TYPE(lri_environment_type), POINTER :: lri_env
251 124119 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
252 : TYPE(mp_para_env_type), POINTER :: para_env
253 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
254 124119 : POINTER :: sab_orb
255 124119 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
256 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
257 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
258 : TYPE(pw_env_type), POINTER :: pw_env
259 : TYPE(pw_poisson_type), POINTER :: poisson_env
260 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
261 372357 : TYPE(pw_r3d_rs_type), DIMENSION(2) :: v_hxc_rspace
262 248238 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace_embed, v_rspace_new, &
263 124119 : v_rspace_new_aux_fit, v_tau_rspace, &
264 124119 : v_tau_rspace_aux_fit
265 : TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs, rho_nlcc, rhoz_cneo_s_rs, v_hartree_rspace, &
266 : v_sccs_rspace, v_sic_rspace, v_spin_ddapc_rest_r, vee, vppl_rspace
267 : TYPE(qs_energy_type), POINTER :: energy
268 124119 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
269 124119 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
270 : TYPE(qs_ks_env_type), POINTER :: ks_env
271 : TYPE(qs_rho_type), POINTER :: rho, rho1, rho_struct, rho_xc
272 : TYPE(section_vals_type), POINTER :: ace_section, &
273 : adiabatic_rescaling_section, &
274 : gauxc_section, hfx_sections, input, &
275 : scf_section, xc_section
276 : TYPE(virial_type), POINTER :: virial
277 :
278 124119 : CALL timeset(routineN, handle)
279 124119 : NULLIFY (admm_env, atomic_kind_set, cell, dft_control, force, logger, mo_derivs, my_rho, &
280 124119 : rho_struct, para_env, pw_env, virial, vppl_rspace, &
281 124119 : ace_section, &
282 124119 : adiabatic_rescaling_section, hfx_sections, input, scf_section, &
283 124119 : xc_section, gauxc_section, matrix_h, matrix_h_im, matrix_s, auxbas_pw_pool, poisson_env, &
284 124119 : v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, matrix_vxc, &
285 124119 : matrix_vxc_kp, matrix_vhxc, &
286 124119 : vee, rho_nlcc, ks_env, ks_matrix, ks_matrix_im, rho, energy, rho_xc, rho_r, rho_ao, &
287 124119 : rho_core, sab_orb, particle_set, qs_kind_set, kpoints)
288 :
289 124119 : CPASSERT(ASSOCIATED(qs_env))
290 :
291 124119 : logger => cp_get_default_logger()
292 124119 : my_print = .TRUE.
293 124119 : IF (PRESENT(print_active)) my_print = print_active
294 124119 : use_gauxc_matrix = .FALSE.
295 124119 : native_gapw_atom_composite_reference = .FALSE.
296 124119 : native_gapw_atom_composite_requested = .FALSE.
297 124119 : native_gapw_composite_direct_ao = .FALSE.
298 124119 : native_gapw_composite_reference = .FALSE.
299 124119 : native_grid_diagnostics = .FALSE.
300 124119 : native_skala_composite_contraction = 0.0_dp
301 124119 : native_skala_composite_rho_contraction = 0.0_dp
302 124119 : native_skala_composite_tau_contraction = 0.0_dp
303 124119 : native_skala_restore_exc = .FALSE.
304 :
305 : CALL get_qs_env(qs_env, &
306 : ks_env=ks_env, &
307 : dft_control=dft_control, &
308 : matrix_h_kp=matrix_h, &
309 : matrix_h_im_kp=matrix_h_im, &
310 : matrix_s_kp=matrix_s, &
311 : matrix_ks_kp=ks_matrix, &
312 : matrix_ks_im_kp=ks_matrix_im, &
313 : matrix_vxc=matrix_vxc, &
314 : matrix_vhxc=matrix_vhxc, &
315 : matrix_vxc_kp=matrix_vxc_kp, &
316 : pw_env=pw_env, &
317 : cell=cell, &
318 : atomic_kind_set=atomic_kind_set, &
319 : para_env=para_env, &
320 : input=input, &
321 : virial=virial, &
322 : v_hartree_rspace=v_hartree_rspace, &
323 : vee=vee, &
324 : rho_nlcc=rho_nlcc, &
325 : rho=rho, &
326 : rho_core=rho_core, &
327 : rho_xc=rho_xc, &
328 : sab_orb=sab_orb, &
329 : energy=energy, &
330 : force=force, &
331 : kpoints=kpoints, &
332 : do_kpoints=do_kpoints, &
333 : particle_set=particle_set, &
334 : qs_kind_set=qs_kind_set, &
335 124119 : natom=natom)
336 :
337 124119 : CALL qs_rho_get(rho, rho_r=rho_r, rho_ao_kp=rho_ao)
338 :
339 124119 : nimages = dft_control%nimages
340 124119 : nspins = dft_control%nspins
341 :
342 : ! remap pointer to allow for non-kpoint external ks matrix
343 124119 : IF (PRESENT(ext_ks_matrix)) ks_matrix(1:nspins, 1:1) => ext_ks_matrix(1:nspins)
344 :
345 124119 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
346 :
347 124119 : adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
348 124119 : CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
349 124119 : just_energy_xc = just_energy
350 124119 : IF (do_adiabatic_rescaling) THEN
351 : !! If we perform adiabatic rescaling, the xc potential has to be scaled by the xc- and
352 : !! HFX-energy. Thus, let us first calculate the energy
353 44 : just_energy_xc = .TRUE.
354 : END IF
355 :
356 124119 : CPASSERT(ASSOCIATED(matrix_h))
357 124119 : CPASSERT(ASSOCIATED(matrix_s))
358 124119 : CPASSERT(ASSOCIATED(rho))
359 124119 : CPASSERT(ASSOCIATED(pw_env))
360 124119 : CPASSERT(SIZE(ks_matrix, 1) > 0)
361 124119 : dokp = (nimages > 1)
362 :
363 : ! Setup the possible usage of DDAPC charges
364 : do_ddapc = dft_control%qs_control%ddapc_restraint .OR. &
365 : qs_env%cp_ddapc_ewald%do_decoupling .OR. &
366 : qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
367 124119 : qs_env%cp_ddapc_ewald%do_solvation
368 :
369 : ! Check if LRIGPW is used
370 124119 : lrigpw = dft_control%qs_control%lrigpw
371 124119 : rigpw = dft_control%qs_control%rigpw
372 124119 : IF (rigpw) THEN
373 26 : CPASSERT(nimages == 1)
374 : END IF
375 26 : IF (lrigpw .AND. rigpw) THEN
376 0 : CPABORT(" LRI and RI are not compatible")
377 : END IF
378 :
379 : ! Check for GAPW method : additional terms for local densities
380 124119 : gapw = dft_control%qs_control%gapw
381 124119 : gapw_xc = dft_control%qs_control%gapw_xc
382 124119 : IF (gapw_xc .AND. gapw) THEN
383 0 : CPABORT(" GAPW and GAPW_XC are not compatible")
384 : END IF
385 124119 : IF ((gapw .AND. lrigpw) .OR. (gapw_xc .AND. lrigpw)) THEN
386 0 : CPABORT(" GAPW/GAPW_XC and LRIGPW are not compatible")
387 : END IF
388 124119 : IF ((gapw .AND. rigpw) .OR. (gapw_xc .AND. rigpw)) THEN
389 0 : CPABORT(" GAPW/GAPW_XC and RIGPW are not compatible")
390 : END IF
391 :
392 124119 : do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
393 124119 : IF (do_ppl) THEN
394 60 : CPASSERT(.NOT. gapw)
395 60 : CALL get_qs_env(qs_env=qs_env, vppl=vppl_rspace)
396 : END IF
397 :
398 124119 : IF (gapw_xc) THEN
399 4336 : CPASSERT(ASSOCIATED(rho_xc))
400 : END IF
401 :
402 : ! gets the tmp grids
403 124119 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
404 :
405 124119 : IF (gapw .AND. (poisson_env%parameters%solver == pw_poisson_implicit)) THEN
406 0 : CPABORT("The implicit Poisson solver cannot be used in conjunction with GAPW.")
407 : END IF
408 :
409 : ! *** Prepare densities for gapw ***
410 124119 : IF (gapw .OR. gapw_xc) THEN
411 26698 : CALL prepare_gapw_den(qs_env, do_rho0=(.NOT. gapw_xc))
412 : END IF
413 :
414 : ! Calculate the Hartree potential
415 124119 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
416 124119 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
417 :
418 124119 : scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
419 : IF (BTEST(cp_print_key_should_output(logger%iter_info, scf_section, &
420 : "PRINT%DETAILED_ENERGY"), &
421 : cp_p_file) .AND. &
422 124119 : (.NOT. gapw) .AND. (.NOT. gapw_xc) .AND. &
423 : (.NOT. (poisson_env%parameters%solver == pw_poisson_implicit))) THEN
424 916 : CALL pw_zero(rho_tot_gspace)
425 916 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density=.TRUE.)
426 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%e_hartree, &
427 916 : v_hartree_gspace)
428 916 : CALL pw_zero(rho_tot_gspace)
429 916 : CALL pw_zero(v_hartree_gspace)
430 : END IF
431 :
432 : ! Get the total density in g-space [ions + electrons]
433 124119 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
434 :
435 124119 : IF (qs_env%scf_control%gce%do_gce .AND. .NOT. dft_control%do_pcc) THEN
436 0 : CPABORT("GCE requires DFT%PLANAR_COUNTER_CHARGE to define the countercharge plane.")
437 : END IF
438 :
439 : ! Add the planar counter charge density
440 124119 : IF (dft_control%do_pcc) THEN
441 114 : CALL planar_counter_charge(rho_tot_gspace, dft_control%pcc_control, auxbas_pw_pool)
442 : END IF
443 :
444 124119 : IF (my_print) THEN
445 124097 : CALL print_densities(qs_env, rho)
446 : END IF
447 :
448 124119 : IF (dft_control%do_sccs) THEN
449 : ! Self-consistent continuum solvation (SCCS) model
450 : NULLIFY (v_sccs_rspace)
451 160 : ALLOCATE (v_sccs_rspace)
452 160 : CALL auxbas_pw_pool%create_pw(v_sccs_rspace)
453 :
454 160 : IF (poisson_env%parameters%solver == pw_poisson_implicit) THEN
455 0 : CPABORT("The implicit Poisson solver cannot be used together with SCCS.")
456 : END IF
457 :
458 160 : IF (use_virial .AND. calculate_forces) THEN
459 : CALL sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs_rspace, &
460 0 : h_stress=h_stress)
461 0 : virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
462 0 : virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
463 : ELSE
464 160 : CALL sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs_rspace)
465 : END IF
466 : ELSE
467 : ! Getting the Hartree energy and Hartree potential. Also getting the stress tensor
468 : ! from the Hartree term if needed. No nuclear force information here
469 123959 : IF (use_virial .AND. calculate_forces) THEN
470 492 : h_stress(:, :) = 0.0_dp
471 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%hartree, &
472 : v_hartree_gspace, h_stress=h_stress, &
473 492 : rho_core=rho_core)
474 6396 : virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
475 6396 : virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
476 : ELSE
477 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%hartree, &
478 123467 : v_hartree_gspace, rho_core=rho_core)
479 : END IF
480 : END IF
481 :
482 124119 : IF (dft_control%do_paep .OR. qs_env%scf_control%gce%do_gce) THEN
483 86 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
484 : CALL planar_averaged_v_hartree_3d(v_hartree_rspace, dft_control, qs_env%scf_control%gce%do_gce, &
485 86 : qs_env%scf_control%gce%ref_esp, para_env)
486 : END IF
487 :
488 : ! In case decouple periodic images and/or apply restraints to charges
489 124119 : IF (do_ddapc) THEN
490 : CALL qs_ks_ddapc(qs_env, auxbas_pw_pool, rho_tot_gspace, v_hartree_gspace, &
491 : v_spin_ddapc_rest_r, energy, calculate_forces, ks_matrix, &
492 1926 : just_energy)
493 : ELSE
494 122193 : dft_control%qs_control%ddapc_explicit_potential = .FALSE.
495 122193 : dft_control%qs_control%ddapc_restraint_is_spin = .FALSE.
496 122193 : IF (.NOT. just_energy) THEN
497 112777 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
498 112777 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
499 : END IF
500 : END IF
501 124119 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
502 :
503 124119 : IF (dft_control%correct_surf_dip) THEN
504 110 : IF (dft_control%surf_dip_correct_switch) THEN
505 110 : CALL calc_dipsurf_potential(qs_env, energy)
506 110 : energy%hartree = energy%hartree + energy%surf_dipole
507 : END IF
508 : END IF
509 :
510 : ! SIC
511 : CALL calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, &
512 124119 : just_energy, calculate_forces, auxbas_pw_pool)
513 :
514 : ! Check if CDFT constraint is needed
515 124119 : CALL qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
516 :
517 : ! Adds the External Potential if requested
518 124119 : IF (dft_control%apply_external_potential) THEN
519 : ! Compute the energy due to the external potential
520 : ee_ener = 0.0_dp
521 728 : DO ispin = 1, nspins
522 728 : ee_ener = ee_ener + pw_integral_ab(rho_r(ispin), vee)
523 : END DO
524 364 : IF (.NOT. just_energy) THEN
525 364 : IF (gapw) THEN
526 : CALL get_qs_env(qs_env=qs_env, &
527 : rho0_s_rs=rho0_s_rs, &
528 42 : rhoz_cneo_s_rs=rhoz_cneo_s_rs)
529 42 : CPASSERT(ASSOCIATED(rho0_s_rs))
530 42 : IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
531 0 : CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs)
532 : END IF
533 42 : ee_ener = ee_ener + pw_integral_ab(rho0_s_rs, vee)
534 42 : IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
535 0 : CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs, -1.0_dp)
536 : END IF
537 : END IF
538 : END IF
539 : ! the sign accounts for the charge of the electrons
540 364 : energy%ee = -ee_ener
541 : END IF
542 :
543 : ! Adds the QM/MM potential
544 124119 : IF (qs_env%qmmm) THEN
545 : CALL qmmm_calculate_energy(qs_env=qs_env, &
546 : rho=rho_r, &
547 : v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, &
548 6318 : qmmm_energy=energy%qmmm_el)
549 6318 : IF (qs_env%qmmm_env_qm%image_charge) THEN
550 : CALL calculate_image_pot(v_hartree_rspace=v_hartree_rspace, &
551 : rho_hartree_gspace=rho_tot_gspace, &
552 : energy=energy, &
553 : qmmm_env=qs_env%qmmm_env_qm, &
554 60 : qs_env=qs_env)
555 60 : IF (.NOT. just_energy) THEN
556 : CALL add_image_pot_to_hartree_pot(v_hartree=v_hartree_rspace, &
557 : v_metal=qs_env%ks_qmmm_env%v_metal_rspace, &
558 60 : qs_env=qs_env)
559 60 : IF (calculate_forces) THEN
560 : CALL integrate_potential_devga_rspace( &
561 : potential=v_hartree_rspace, coeff=qs_env%image_coeff, &
562 : forces=qs_env%qmmm_env_qm%image_charge_pot%image_forcesMM, &
563 20 : qmmm_env=qs_env%qmmm_env_qm, qs_env=qs_env)
564 : END IF
565 : END IF
566 60 : CALL qs_env%ks_qmmm_env%v_metal_rspace%release()
567 60 : DEALLOCATE (qs_env%ks_qmmm_env%v_metal_rspace)
568 : END IF
569 6318 : IF (.NOT. just_energy) THEN
570 : CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, &
571 6228 : v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, scale=1.0_dp)
572 : END IF
573 : END IF
574 124119 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
575 :
576 : ! SMEAGOL interface
577 124119 : IF (dft_control%smeagol_control%smeagol_enabled .AND. &
578 : dft_control%smeagol_control%run_type == smeagol_runtype_emtransport) THEN
579 0 : CPASSERT(ASSOCIATED(dft_control%smeagol_control%aux))
580 : CALL smeagol_shift_v_hartree(v_hartree_rspace, cell, &
581 : dft_control%smeagol_control%aux%HartreeLeadsLeft, &
582 : dft_control%smeagol_control%aux%HartreeLeadsRight, &
583 : dft_control%smeagol_control%aux%HartreeLeadsBottom, &
584 : dft_control%smeagol_control%aux%VBias, &
585 : dft_control%smeagol_control%aux%minL, &
586 : dft_control%smeagol_control%aux%maxR, &
587 : dft_control%smeagol_control%aux%isexplicit_maxR, &
588 0 : dft_control%smeagol_control%aux%isexplicit_HartreeLeadsBottom)
589 : END IF
590 :
591 : ! calculate the density matrix for the fitted mo_coeffs
592 124119 : IF (dft_control%do_admm) THEN
593 13290 : IF (PRESENT(ext_xc_section)) THEN
594 0 : CALL hfx_admm_init(qs_env, calculate_forces, ext_xc_section)
595 : ELSE
596 13290 : CALL hfx_admm_init(qs_env, calculate_forces)
597 : END IF
598 :
599 13290 : IF (dft_control%do_admm_mo) THEN
600 13076 : IF (qs_env%run_rtp) THEN
601 92 : CALL rtp_admm_calc_rho_aux(qs_env)
602 : ELSE
603 12984 : IF (dokp) THEN
604 156 : CALL admm_mo_calc_rho_aux_kp(qs_env)
605 : ELSE
606 12828 : CALL admm_mo_calc_rho_aux(qs_env)
607 : END IF
608 : END IF
609 214 : ELSE IF (dft_control%do_admm_dm) THEN
610 214 : CALL admm_dm_calc_rho_aux(qs_env)
611 : END IF
612 : END IF
613 :
614 : ! only activate stress calculation if
615 124119 : IF (use_virial .AND. calculate_forces) virial%pv_calculate = .TRUE.
616 :
617 : ! *** calculate the xc potential on the pw density ***
618 : ! *** associates v_rspace_new if the xc potential needs to be computed.
619 : ! If we do wavefunction fitting, we need the vxc_potential in the auxiliary basis set
620 124119 : IF (dft_control%do_admm) THEN
621 13290 : CALL get_qs_env(qs_env, admm_env=admm_env)
622 13290 : xc_section => admm_env%xc_section_aux
623 13290 : CALL get_admm_env(admm_env, rho_aux_fit=rho_struct)
624 :
625 : ! here we ignore a possible vdW section in admm_env%xc_section_aux
626 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
627 : vxc_rho=v_rspace_new_aux_fit, vxc_tau=v_tau_rspace_aux_fit, exc=energy%exc_aux_fit, &
628 13290 : just_energy=just_energy_xc)
629 :
630 13290 : IF (admm_env%do_gapw) THEN
631 : !compute the potential due to atomic densities
632 : CALL calculate_vxc_atom(qs_env, energy_only=just_energy_xc, exc1=energy%exc1_aux_fit, &
633 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
634 : xc_section_external=xc_section, &
635 : rho_atom_set_external=admm_env%admm_gapw_env%local_rho_set%rho_atom_set, &
636 4674 : calculate_forces=calculate_forces)
637 :
638 : END IF
639 :
640 13290 : NULLIFY (rho_struct)
641 :
642 13290 : IF (use_virial .AND. calculate_forces) THEN
643 20 : vscale = 1.0_dp
644 : !Note: ADMMS and ADMMP stress tensor only for closed-shell calculations
645 20 : IF (admm_env%do_admms) vscale = admm_env%gsi(1)**(2.0_dp/3.0_dp)
646 20 : IF (admm_env%do_admmp) vscale = admm_env%gsi(1)**2
647 260 : virial%pv_exc = virial%pv_exc - vscale*virial%pv_xc
648 260 : virial%pv_virial = virial%pv_virial - vscale*virial%pv_xc
649 : ! virial%pv_xc will be zeroed in the xc routines
650 : END IF
651 13290 : xc_section => admm_env%xc_section_primary
652 : ELSE
653 110829 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
654 : ! build ks matrix with an xc section potentially different from the one defined in input
655 110829 : IF (PRESENT(ext_xc_section)) xc_section => ext_xc_section
656 : END IF
657 :
658 124119 : IF (gapw_xc) THEN
659 4336 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
660 : ELSE
661 119783 : CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
662 : END IF
663 :
664 : ! zmp
665 124119 : IF (dft_control%apply_external_density .OR. dft_control%apply_external_vxc) THEN
666 0 : energy%exc = 0.0_dp
667 0 : CALL calculate_zmp_potential(qs_env, v_rspace_new, rho, exc=energy%exc)
668 : ELSE
669 : ! Embedding potential (runs regardless of XC method)
670 124119 : IF (dft_control%apply_embed_pot) THEN
671 868 : NULLIFY (v_rspace_embed)
672 868 : energy%embed_corr = 0.0_dp
673 : CALL get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, &
674 868 : energy%embed_corr, just_energy)
675 : END IF
676 :
677 : ! Everything else, either via GauXC or manual XC computation
678 124119 : native_skala_evaluator = xc_section_uses_native_skala_grid(xc_section)
679 124119 : paw_composite_representation = .FALSE.
680 124119 : IF (dft_control%use_gauxc .AND. (gapw .OR. gapw_xc)) THEN
681 : paw_composite_representation = &
682 : skala_gapw_representation(xc_section) == skala_gapw_paw_one_center .AND. &
683 286 : gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
684 : native_skala_evaluator = native_skala_evaluator .OR. &
685 286 : paw_composite_representation
686 : END IF
687 124119 : IF (dft_control%use_gauxc) THEN
688 804 : IF (native_skala_evaluator) THEN
689 326 : CALL ensure_native_skala_grid_scope(xc_section)
690 : native_gapw_composite_reference = &
691 326 : native_skala_gapw_composite_reference(xc_section) .AND. (gapw .OR. gapw_xc)
692 : native_gapw_atom_composite_requested = &
693 : native_skala_gapw_atom_composite_requested(xc_section) .AND. &
694 326 : (gapw .OR. gapw_xc)
695 : native_gapw_atom_composite_reference = &
696 326 : native_gapw_atom_composite_requested .OR. paw_composite_representation
697 458 : IF (native_gapw_atom_composite_requested .AND. ANY(cell%perd /= 0)) THEN
698 : CALL cp_abort(__LOCATION__, &
699 0 : "The atom-centered GAPW composite reference is molecular only.")
700 : END IF
701 326 : IF (paw_composite_representation) THEN
702 208 : IF (ANY(cell%perd /= 0)) THEN
703 96 : native_gapw_composite_reference = .TRUE.
704 96 : native_gapw_atom_composite_reference = .FALSE.
705 : ELSE
706 28 : native_gapw_composite_reference = .FALSE.
707 28 : native_gapw_atom_composite_reference = .TRUE.
708 : END IF
709 : END IF
710 326 : IF (native_gapw_composite_reference .AND. native_gapw_atom_composite_reference) THEN
711 : CALL cp_abort(__LOCATION__, &
712 0 : "Select only one native-grid GAPW composite reference route.")
713 : END IF
714 : native_gapw_composite_direct_ao = native_gapw_composite_reference .AND. &
715 326 : native_skala_gapw_composite_direct_ao(xc_section)
716 326 : gauxc_section => get_gauxc_section(xc_section)
717 326 : IF (ASSOCIATED(gauxc_section)) THEN
718 : CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_DIAGNOSTICS", &
719 326 : l_val=native_grid_diagnostics)
720 : END IF
721 326 : IF ((.NOT. do_kpoints) .AND. nimages /= 1) THEN
722 : CALL cp_abort(__LOCATION__, &
723 : "Native SKALA grid evaluation supports multiple images only "// &
724 0 : "for k-point calculations.")
725 : END IF
726 326 : IF (do_kpoints) THEN
727 48 : CPASSERT(ASSOCIATED(kpoints))
728 48 : gauxc_section => get_gauxc_section(xc_section)
729 48 : CPASSERT(ASSOCIATED(gauxc_section))
730 48 : CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_USE_CUDA", l_val=native_grid_use_cuda)
731 48 : IF (.NOT. native_grid_use_cuda) THEN
732 48 : IF (para_env%mepos == 0 .AND. .NOT. native_grid_cpu_kpoints_warned) THEN
733 : CALL cp_warn(__LOCATION__, &
734 : "Native SKALA grid evaluation with k-points is using the CPU TorchScript "// &
735 : "path. Use a mutually compatible BLAS, ScaLAPACK, OpenMP, and LibTorch "// &
736 12 : "runtime stack, or select NATIVE_GRID_USE_CUDA T.")
737 12 : native_grid_cpu_kpoints_warned = .TRUE.
738 : END IF
739 : END IF
740 : END IF
741 326 : IF (dft_control%roks) THEN
742 0 : CPABORT("Native SKALA grid evaluation does not support ROKS.")
743 : END IF
744 326 : IF (dft_control%do_admm) THEN
745 0 : CPABORT("Native SKALA grid evaluation does not support ADMM.")
746 : END IF
747 : ! Force and stress rebuilds re-enter this path for derivatives and VXC only. The
748 : ! autograd virial already contains the rebuilt XC volume term when qs_vxc_create
749 : ! returns, so the externally reported energy must be restored to its SCF value.
750 326 : native_skala_restore_exc = calculate_forces
751 326 : IF (native_skala_restore_exc) THEN
752 68 : native_skala_exc_scf = energy%exc
753 68 : native_skala_total_scf = energy%total
754 : END IF
755 326 : IF (calculate_forces) THEN
756 204 : ALLOCATE (native_skala_atom_force(3, natom))
757 68 : native_skala_atom_force = 0.0_dp
758 68 : IF (native_gapw_atom_composite_reference) THEN
759 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
760 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
761 : edisp=edisp, dispersion_env=qs_env%dispersion_env, &
762 : just_energy=just_energy_xc, qs_env_external=qs_env, &
763 : native_gapw_composite_override=native_gapw_composite_reference, &
764 4 : native_skala_defer_to_atom_composite=.TRUE.)
765 : ELSE
766 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
767 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
768 : edisp=edisp, dispersion_env=qs_env%dispersion_env, &
769 : just_energy=just_energy_xc, &
770 : native_skala_atom_force=native_skala_atom_force, &
771 : qs_env_external=qs_env, &
772 64 : native_gapw_composite_override=native_gapw_composite_reference)
773 : END IF
774 : IF (native_gapw_composite_reference .AND. .NOT. native_gapw_composite_direct_ao .AND. &
775 68 : ASSOCIATED(v_rspace_new) .AND. ASSOCIATED(v_tau_rspace)) THEN
776 20 : IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
777 1 : output_unit = cp_logger_get_default_io_unit()
778 1 : IF (output_unit > 0) THEN
779 3 : DO iatom = 1, natom
780 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,I0,3(1X,ES20.12))") &
781 2 : "SKALA_GPW| Native regular-grid atom force", iatom, &
782 5 : native_skala_atom_force(:, iatom)
783 : END DO
784 : END IF
785 : END IF
786 20 : IF (use_virial) THEN
787 20 : native_skala_composite_virial = 0.0_dp
788 : CALL integrate_gapw_composite_vxc_gspace( &
789 : qs_env, auxbas_pw_pool, v_rspace_new, v_tau_rspace, &
790 : atom_force=native_skala_atom_force, &
791 : strain_virial=native_skala_composite_virial, &
792 : one_center_contraction=native_skala_composite_contraction, &
793 : one_center_rho_contraction=native_skala_composite_rho_contraction, &
794 20 : one_center_tau_contraction=native_skala_composite_tau_contraction)
795 20 : IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
796 1 : output_unit = cp_logger_get_default_io_unit()
797 1 : IF (output_unit > 0) THEN
798 : WRITE (UNIT=output_unit, FMT="(T2,A)") &
799 1 : "SKALA_GPW| Composite-direction one-center virial"
800 4 : DO iatom = 1, 3
801 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,3ES20.10)") &
802 4 : "SKALA_GPW|", native_skala_composite_virial(iatom, :)
803 : END DO
804 : END IF
805 : END IF
806 260 : virial%pv_xc = virial%pv_xc + native_skala_composite_virial
807 : ELSE
808 : CALL integrate_gapw_composite_vxc_gspace( &
809 : qs_env, auxbas_pw_pool, v_rspace_new, v_tau_rspace, &
810 : atom_force=native_skala_atom_force, &
811 : one_center_contraction=native_skala_composite_contraction, &
812 : one_center_rho_contraction=native_skala_composite_rho_contraction, &
813 0 : one_center_tau_contraction=native_skala_composite_tau_contraction)
814 : END IF
815 : END IF
816 68 : IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
817 4 : output_unit = cp_logger_get_default_io_unit()
818 4 : IF (output_unit > 0) THEN
819 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
820 4 : "SKALA_GPW| Composite-direction one-center VXC contraction", &
821 8 : native_skala_composite_contraction
822 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
823 4 : "SKALA_GPW| Composite-direction one-center rho contraction", &
824 8 : native_skala_composite_rho_contraction
825 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
826 4 : "SKALA_GPW| Composite-direction one-center tau contraction", &
827 8 : native_skala_composite_tau_contraction
828 12 : DO iatom = 1, natom
829 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,I0,3(1X,ES20.12))") &
830 12 : "SKALA_GPW| Native atom force", iatom, native_skala_atom_force(:, iatom)
831 : END DO
832 : END IF
833 : END IF
834 68 : CPASSERT(ASSOCIATED(force))
835 68 : CPASSERT(ASSOCIATED(atomic_kind_set))
836 68 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
837 204 : DO iatom = 1, natom
838 136 : ikind = kind_of(iatom)
839 136 : atom_a = atom_of_kind(iatom)
840 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + &
841 612 : native_skala_atom_force(:, iatom)
842 : END DO
843 68 : DEALLOCATE (atom_of_kind, kind_of, native_skala_atom_force)
844 : ELSE
845 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
846 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
847 : edisp=edisp, dispersion_env=qs_env%dispersion_env, &
848 : just_energy=just_energy_xc, qs_env_external=qs_env, &
849 : native_gapw_composite_override=native_gapw_composite_reference, &
850 : native_skala_defer_to_atom_composite= &
851 258 : native_gapw_atom_composite_reference)
852 : IF (native_gapw_composite_reference .AND. .NOT. native_gapw_composite_direct_ao .AND. &
853 258 : ASSOCIATED(v_rspace_new) .AND. ASSOCIATED(v_tau_rspace)) THEN
854 : CALL integrate_gapw_composite_vxc_gspace( &
855 : qs_env, auxbas_pw_pool, v_rspace_new, v_tau_rspace, &
856 : one_center_contraction=native_skala_composite_contraction, &
857 : one_center_rho_contraction=native_skala_composite_rho_contraction, &
858 76 : one_center_tau_contraction=native_skala_composite_tau_contraction)
859 76 : IF (native_grid_diagnostics .AND. para_env%mepos == 0) THEN
860 1 : output_unit = cp_logger_get_default_io_unit()
861 1 : IF (output_unit > 0) THEN
862 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
863 1 : "SKALA_GPW| Composite-direction one-center VXC contraction", &
864 2 : native_skala_composite_contraction
865 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
866 1 : "SKALA_GPW| Composite-direction one-center rho contraction", &
867 2 : native_skala_composite_rho_contraction
868 : WRITE (UNIT=output_unit, FMT="(T2,A,1X,ES20.10)") &
869 1 : "SKALA_GPW| Composite-direction one-center tau contraction", &
870 2 : native_skala_composite_tau_contraction
871 : END IF
872 : END IF
873 : END IF
874 : END IF
875 326 : IF (native_skala_restore_exc) energy%exc = native_skala_exc_scf
876 326 : IF ((gapw .OR. gapw_xc) .AND. .NOT. native_gapw_composite_reference) THEN
877 80 : IF (native_gapw_atom_composite_reference) THEN
878 28 : CPASSERT(ASSOCIATED(v_rspace_new))
879 28 : CPASSERT(ASSOCIATED(v_tau_rspace))
880 : CALL calculate_vxc_atom( &
881 : qs_env, just_energy_xc, energy%exc1, xc_section_external=xc_section, &
882 : calculate_forces=calculate_forces, composite_vxc_rho=v_rspace_new, &
883 28 : composite_vxc_tau=v_tau_rspace)
884 28 : energy%exc = 0.0_dp
885 : ELSE
886 : CALL calculate_vxc_atom(qs_env, just_energy_xc, energy%exc1, &
887 : xc_section_external=xc_section, &
888 52 : calculate_forces=calculate_forces)
889 : END IF
890 : END IF
891 326 : IF (edisp /= 0.0_dp) energy%dispersion = edisp
892 326 : IF (qs_env%requires_matrix_vxc .AND. ASSOCIATED(v_rspace_new)) THEN
893 0 : IF (do_kpoints) THEN
894 : CALL compute_matrix_vxc_kp(qs_env=qs_env, v_rspace=v_rspace_new, &
895 : matrix_vxc_kp=matrix_vxc_kp, &
896 0 : gapw_full_basis=native_gapw_composite_direct_ao)
897 0 : CALL set_ks_env(ks_env, matrix_vxc_kp=matrix_vxc_kp)
898 : ELSE
899 : CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace_new, &
900 : matrix_vxc=matrix_vxc, &
901 0 : gapw_full_basis=native_gapw_composite_direct_ao)
902 0 : CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
903 : END IF
904 : END IF
905 : ELSE
906 478 : use_gauxc_matrix = .TRUE.
907 478 : CALL apply_gauxc(qs_env, xc_section, calculate_forces)
908 478 : IF (gapw_xc .OR. (gapw .AND. gauxc_gapw_has_paw_pseudopotentials(qs_kind_set))) THEN
909 : CALL calculate_vxc_atom(qs_env, just_energy_xc, energy%exc1, &
910 : xc_section_external=xc_section, &
911 2 : calculate_forces=calculate_forces)
912 : END IF
913 : END IF
914 : ELSE
915 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
916 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
917 : edisp=edisp, dispersion_env=qs_env%dispersion_env, &
918 123315 : just_energy=just_energy_xc)
919 123315 : IF (edisp /= 0.0_dp) energy%dispersion = edisp
920 123315 : IF (qs_env%requires_matrix_vxc .AND. ASSOCIATED(v_rspace_new)) THEN
921 2 : CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace_new, matrix_vxc=matrix_vxc)
922 2 : CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
923 : END IF
924 :
925 123315 : IF (gapw .OR. gapw_xc) THEN
926 : CALL calculate_vxc_atom(qs_env, just_energy_xc, energy%exc1, &
927 : xc_section_external=xc_section, &
928 26412 : calculate_forces=calculate_forces)
929 : END IF
930 : END IF
931 : END IF
932 :
933 26848 : IF (gapw .AND. dft_control%qs_control%cdft) THEN
934 1260 : ALLOCATE (cdft_gapw_values(SIZE(cdft_control%group)))
935 420 : IF (cdft_control%atomic_charges) THEN
936 : ALLOCATE (cdft_charge_correction(natom, nspins), &
937 532 : cdft_electronic_charge(cdft_control%natoms, nspins))
938 : CALL gapw_cdft_one_center(qs_env, just_energy, calculate_forces, &
939 76 : cdft_gapw_values, cdft_charge_correction)
940 228 : DO ispin = 1, nspins
941 624 : DO iatom = 1, cdft_control%natoms
942 : cdft_electronic_charge(iatom, ispin) = &
943 548 : pw_integral_ab(cdft_control%charge(iatom), rho_r(ispin), local_only=.TRUE.)
944 : END DO
945 : END DO
946 76 : CALL para_env%sum(cdft_electronic_charge)
947 228 : DO ispin = 1, nspins
948 624 : DO iatom = 1, cdft_control%natoms
949 : cdft_electronic_charge(iatom, ispin) = &
950 : cdft_electronic_charge(iatom, ispin) + &
951 548 : cdft_charge_correction(cdft_control%atoms(iatom), ispin)
952 : END DO
953 : END DO
954 : ELSE
955 344 : CALL gapw_cdft_one_center(qs_env, just_energy, calculate_forces, cdft_gapw_values)
956 : END IF
957 912 : cdft_control%value = cdft_control%value + cdft_gapw_values
958 912 : energy%cdft = energy%cdft + DOT_PRODUCT(cdft_control%strength, cdft_gapw_values)
959 420 : IF (cdft_control%atomic_charges) THEN
960 76 : CALL cdft_constraint_print(qs_env, cdft_electronic_charge)
961 76 : DEALLOCATE (cdft_charge_correction, cdft_electronic_charge)
962 : ELSE
963 344 : CALL cdft_constraint_print(qs_env, cdft_electronic_charge)
964 : END IF
965 420 : DEALLOCATE (cdft_gapw_values)
966 : END IF
967 :
968 124119 : IF (dft_control%mtlr_dft_with_perturbation .AND. .NOT. just_energy) THEN
969 420 : IF (do_adiabatic_rescaling .OR. use_gauxc_matrix) THEN
970 0 : CPABORT("MTLR requires an XC potential on the real-space grid.")
971 : END IF
972 420 : IF (.NOT. ASSOCIATED(matrix_vhxc)) THEN
973 20 : ALLOCATE (matrix_vhxc(nspins))
974 12 : DO ispin = 1, nspins
975 8 : NULLIFY (matrix_vhxc(ispin)%matrix)
976 8 : ALLOCATE (matrix_vhxc(ispin)%matrix)
977 : CALL dbcsr_create(matrix_vhxc(ispin)%matrix, &
978 8 : template=ks_matrix(1, 1)%matrix)
979 8 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_vhxc(ispin)%matrix, sab_orb)
980 12 : CALL dbcsr_set(matrix_vhxc(ispin)%matrix, 0.0_dp)
981 : END DO
982 4 : CALL set_ks_env(ks_env, matrix_vhxc=matrix_vhxc)
983 : ELSE
984 416 : CPASSERT(SIZE(matrix_vhxc) == nspins)
985 1248 : DO ispin = 1, nspins
986 1248 : CPASSERT(ASSOCIATED(matrix_vhxc(ispin)%matrix))
987 : END DO
988 : END IF
989 1260 : DO ispin = 1, nspins
990 840 : CALL auxbas_pw_pool%create_pw(v_hxc_rspace(ispin))
991 840 : CALL pw_zero(v_hxc_rspace(ispin))
992 840 : CALL pw_axpy(v_hartree_rspace, v_hxc_rspace(ispin))
993 840 : IF (ASSOCIATED(v_rspace_new)) THEN
994 : CALL pw_axpy(v_rspace_new(ispin), &
995 : v_hxc_rspace(ispin), &
996 840 : v_rspace_new(ispin)%pw_grid%dvol)
997 : END IF
998 840 : CALL dbcsr_set(matrix_vhxc(ispin)%matrix, 0.0_dp)
999 : CALL integrate_v_rspace(v_hxc_rspace(ispin), &
1000 : hmat=matrix_vhxc(ispin), &
1001 : qs_env=qs_env, &
1002 840 : calculate_forces=.FALSE.)
1003 1260 : CALL auxbas_pw_pool%give_back_pw(v_hxc_rspace(ispin))
1004 : END DO
1005 : END IF
1006 :
1007 : ! set hartree and xc potentials for use in Harris method
1008 124119 : IF (qs_env%harris_method) THEN
1009 80 : CALL get_qs_env(qs_env, harris_env=harris_env)
1010 80 : CALL harris_set_potentials(harris_env, v_hartree_rspace, v_rspace_new)
1011 : END IF
1012 :
1013 124119 : NULLIFY (rho_struct)
1014 124119 : IF (use_virial .AND. calculate_forces) THEN
1015 6396 : virial%pv_exc = virial%pv_exc - virial%pv_xc
1016 6396 : virial%pv_virial = virial%pv_virial - virial%pv_xc
1017 : END IF
1018 :
1019 : ! *** Add Hartree-Fock contribution if required ***
1020 124119 : hfx_sections => section_vals_get_subs_vals(xc_section, "HF")
1021 124119 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
1022 :
1023 124119 : ace_active = .FALSE.
1024 124119 : ace_rebuild_frequency = 1
1025 :
1026 124119 : IF (do_hfx) THEN
1027 29634 : ace_section => section_vals_get_subs_vals(hfx_sections, "ACE")
1028 29634 : IF (ASSOCIATED(ace_section)) THEN
1029 29634 : CALL section_vals_val_get(ace_section, "ACTIVE", l_val=ace_active)
1030 29634 : CALL section_vals_val_get(ace_section, "REBUILD_FREQUENCY", i_val=ace_rebuild_frequency)
1031 : END IF
1032 : END IF
1033 :
1034 124119 : IF (do_hfx) THEN
1035 29634 : IF (dokp) THEN
1036 274 : IF (ace_active) THEN
1037 0 : CPABORT("ACE-HFX for k-points is not implemented yet")
1038 : ELSE
1039 274 : CALL hfx_ks_matrix_kp(qs_env, ks_matrix, energy, calculate_forces)
1040 : END IF
1041 :
1042 : ELSE
1043 : ! ext_xc_section may contain a hfx section
1044 29360 : IF (ace_active) THEN
1045 : CALL hfx_ace_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, &
1046 : just_energy, v_rspace_new, v_tau_rspace, &
1047 48 : ace_rebuild_frequency, ext_xc_section=xc_section)
1048 : ELSE
1049 : CALL hfx_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, &
1050 29312 : just_energy, v_rspace_new, v_tau_rspace, ext_xc_section=xc_section)
1051 : END IF
1052 : END IF
1053 : END IF !do_hfx
1054 :
1055 124119 : IF (do_ppl .AND. calculate_forces) THEN
1056 12 : CPASSERT(.NOT. gapw)
1057 26 : DO ispin = 1, nspins
1058 26 : CALL integrate_ppl_rspace(rho_r(ispin), qs_env)
1059 : END DO
1060 : END IF
1061 :
1062 : ! The atom-composite evaluator differentiates its atom-centered NLCC field directly.
1063 124119 : IF (ASSOCIATED(rho_nlcc) .AND. calculate_forces .AND. &
1064 : .NOT. native_gapw_atom_composite_reference) THEN
1065 80 : DO ispin = 1, nspins
1066 40 : CALL integrate_rho_nlcc(v_rspace_new(ispin), qs_env)
1067 80 : IF (dft_control%do_admm) CALL integrate_rho_nlcc(v_rspace_new_aux_fit(ispin), qs_env)
1068 : END DO
1069 : END IF
1070 :
1071 : ! calculate KG correction
1072 124119 : IF (dft_control%qs_control%do_kg .AND. just_energy) THEN
1073 :
1074 12 : CPASSERT(nimages == 1)
1075 12 : ksmat => ks_matrix(:, 1)
1076 12 : CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces, do_kernel=.FALSE.)
1077 :
1078 : ! subtract kg corr from the total energy
1079 12 : energy%exc = energy%exc - ekin_mol
1080 :
1081 : END IF
1082 :
1083 : ! *** Single atom contributions ***
1084 124119 : IF (.NOT. just_energy) THEN
1085 114313 : IF (calculate_forces) THEN
1086 : ! Getting nuclear force contribution from the core charge density
1087 6013 : IF ((poisson_env%parameters%solver == pw_poisson_implicit) .AND. &
1088 : (poisson_env%parameters%dielectric_params%dielec_core_correction)) THEN
1089 28 : BLOCK
1090 : TYPE(pw_r3d_rs_type) :: v_minus_veps
1091 28 : CALL auxbas_pw_pool%create_pw(v_minus_veps)
1092 28 : CALL pw_copy(v_hartree_rspace, v_minus_veps)
1093 28 : CALL pw_axpy(poisson_env%implicit_env%v_eps, v_minus_veps, -v_hartree_rspace%pw_grid%dvol)
1094 28 : CALL integrate_v_core_rspace(v_minus_veps, qs_env)
1095 28 : CALL auxbas_pw_pool%give_back_pw(v_minus_veps)
1096 : END BLOCK
1097 : ELSE
1098 5985 : CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
1099 : END IF
1100 : END IF
1101 :
1102 114313 : IF (.NOT. do_hfx) THEN
1103 : ! Initialize the Kohn-Sham matrix with the core Hamiltonian matrix
1104 : ! (sets ks sparsity equal to matrix_h sparsity)
1105 190129 : DO ispin = 1, nspins
1106 700039 : DO img = 1, nimages
1107 509910 : CALL dbcsr_get_info(ks_matrix(ispin, img)%matrix, name=name) ! keep the name
1108 612984 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix, name=name)
1109 : END DO
1110 : END DO
1111 : ! imaginary part if required
1112 87055 : IF (qs_env%run_rtp) THEN
1113 2036 : IF (dft_control%rtp_control%velocity_gauge) THEN
1114 150 : CPASSERT(ASSOCIATED(matrix_h_im))
1115 150 : CPASSERT(ASSOCIATED(ks_matrix_im))
1116 300 : DO ispin = 1, nspins
1117 450 : DO img = 1, nimages
1118 150 : CALL dbcsr_get_info(ks_matrix_im(ispin, img)%matrix, name=name) ! keep the name
1119 300 : CALL dbcsr_copy(ks_matrix_im(ispin, img)%matrix, matrix_h_im(1, img)%matrix, name=name)
1120 : END DO
1121 : END DO
1122 : END IF
1123 : END IF
1124 : END IF
1125 :
1126 114313 : IF (use_virial .AND. calculate_forces) THEN
1127 6396 : pv_loc = virial%pv_virial
1128 : END IF
1129 : ! sum up potentials and integrate
1130 : ! Pointing my_rho to the density matrix rho_ao
1131 114313 : my_rho => rho_ao
1132 :
1133 : CALL sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, &
1134 : v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, &
1135 : v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, &
1136 114313 : cdft_control, calculate_forces)
1137 :
1138 114313 : IF (use_gauxc_matrix) THEN
1139 478 : IF (dokp) THEN
1140 0 : CALL get_qs_env(qs_env=qs_env, matrix_vxc_kp=matrix_vxc_kp)
1141 0 : CPASSERT(ASSOCIATED(matrix_vxc_kp))
1142 0 : DO ispin = 1, nspins
1143 0 : DO img = 1, nimages
1144 : CALL dbcsr_add(ks_matrix(ispin, img)%matrix, matrix_vxc_kp(ispin, img)%matrix, &
1145 0 : 1.0_dp, 1.0_dp)
1146 : END DO
1147 : END DO
1148 : ELSE
1149 478 : CALL get_qs_env(qs_env=qs_env, matrix_vxc=matrix_vxc)
1150 478 : CPASSERT(ASSOCIATED(matrix_vxc))
1151 478 : CPASSERT(nimages == 1)
1152 978 : DO ispin = 1, nspins
1153 978 : CALL dbcsr_add(ks_matrix(ispin, 1)%matrix, matrix_vxc(ispin)%matrix, 1.0_dp, 1.0_dp)
1154 : END DO
1155 : END IF
1156 : END IF
1157 :
1158 114313 : IF (gapw .OR. gapw_xc) THEN
1159 25808 : IF (calculate_forces) THEN
1160 774 : IF (gapw_xc) THEN
1161 124 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
1162 : ELSE
1163 650 : CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
1164 : END IF
1165 774 : NULLIFY (rho1)
1166 774 : IF (dft_control%use_gauxc .AND. (gapw .OR. gapw_xc) .AND. &
1167 : .NOT. native_skala_evaluator) THEN
1168 2 : gauxc_model_none = .FALSE.
1169 2 : gauxc_section => get_gauxc_section(xc_section)
1170 2 : IF (ASSOCIATED(gauxc_section)) THEN
1171 2 : CALL section_vals_val_get(gauxc_section, "MODEL", c_val=gauxc_model_name)
1172 2 : gauxc_model_name = ADJUSTL(gauxc_model_name)
1173 2 : CALL uppercase(gauxc_model_name)
1174 : gauxc_model_none = (TRIM(gauxc_model_name) == "" .OR. &
1175 2 : TRIM(gauxc_model_name) == "NONE")
1176 : END IF
1177 2 : IF (gauxc_model_none .AND. &
1178 : (gapw_xc .OR. gauxc_gapw_has_paw_pseudopotentials(qs_kind_set))) THEN
1179 0 : CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section)
1180 : END IF
1181 : ELSE
1182 772 : CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section)
1183 : END IF
1184 : !
1185 774 : IF (dft_control%do_admm) THEN
1186 90 : CALL get_qs_env(qs_env, admm_env=admm_env)
1187 90 : xc_section => admm_env%xc_section_aux
1188 90 : CALL get_admm_env(admm_env, rho_aux_fit=rho_struct)
1189 90 : vscale = 1.0_dp
1190 90 : IF (admm_env%do_admmp) THEN
1191 8 : vscale = admm_env%gsi(1)**2
1192 82 : ELSE IF (admm_env%do_admms) THEN
1193 6 : vscale = admm_env%gsi(1)**(2.0_dp/3.0_dp)
1194 : END IF
1195 90 : CALL accint_weight_force(qs_env, rho_struct, rho1, 0, xc_section, force_scale=vscale)
1196 : END IF
1197 : END IF
1198 : END IF
1199 :
1200 114313 : IF (use_virial .AND. calculate_forces) THEN
1201 6396 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1202 : END IF
1203 114313 : IF (dft_control%qs_control%do_kg) THEN
1204 978 : CPASSERT(nimages == 1)
1205 978 : ksmat => ks_matrix(:, 1)
1206 :
1207 978 : IF (use_virial .AND. calculate_forces) THEN
1208 208 : pv_loc = virial%pv_virial
1209 : END IF
1210 :
1211 978 : CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces, do_kernel=.FALSE.)
1212 : ! subtract kg corr from the total energy
1213 978 : energy%exc = energy%exc - ekin_mol
1214 :
1215 : ! virial corrections
1216 978 : IF (use_virial .AND. calculate_forces) THEN
1217 :
1218 : ! Integral contribution
1219 208 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
1220 :
1221 : ! GGA contribution
1222 208 : virial%pv_exc = virial%pv_exc + virial%pv_xc
1223 208 : virial%pv_virial = virial%pv_virial + virial%pv_xc
1224 208 : virial%pv_xc = 0.0_dp
1225 : END IF
1226 : END IF
1227 :
1228 : ELSE
1229 : IF (do_hfx) THEN
1230 : IF (.FALSE.) THEN
1231 : CPWARN("KS matrix no longer correct. Check possible problems with property calculations!")
1232 : END IF
1233 : END IF
1234 : END IF ! .NOT. just energy
1235 124119 : IF (dft_control%qs_control%ddapc_explicit_potential) THEN
1236 164 : CALL auxbas_pw_pool%give_back_pw(v_spin_ddapc_rest_r)
1237 164 : DEALLOCATE (v_spin_ddapc_rest_r)
1238 : END IF
1239 :
1240 124119 : IF (calculate_forces .AND. dft_control%qs_control%cdft) THEN
1241 126 : IF (.NOT. cdft_control%transfer_pot) THEN
1242 230 : DO iatom = 1, SIZE(cdft_control%group)
1243 124 : CALL auxbas_pw_pool%give_back_pw(cdft_control%group(iatom)%weight)
1244 230 : DEALLOCATE (cdft_control%group(iatom)%weight)
1245 : END DO
1246 106 : IF (cdft_control%atomic_charges) THEN
1247 78 : DO iatom = 1, cdft_control%natoms
1248 78 : CALL auxbas_pw_pool%give_back_pw(cdft_control%charge(iatom))
1249 : END DO
1250 26 : DEALLOCATE (cdft_control%charge)
1251 : END IF
1252 106 : IF (cdft_control%type == outer_scf_becke_constraint .AND. &
1253 : cdft_control%becke_control%cavity_confine) THEN
1254 88 : IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN
1255 64 : CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
1256 : ELSE
1257 24 : DEALLOCATE (cdft_control%becke_control%cavity_mat)
1258 : END IF
1259 18 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
1260 6 : IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) THEN
1261 0 : CALL auxbas_pw_pool%give_back_pw(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)
1262 : END IF
1263 : END IF
1264 106 : IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment)
1265 106 : cdft_control%save_pot = .FALSE.
1266 106 : cdft_control%need_pot = .TRUE.
1267 106 : cdft_control%external_control = .FALSE.
1268 : END IF
1269 : END IF
1270 :
1271 124119 : IF (dft_control%do_sccs) THEN
1272 160 : CALL auxbas_pw_pool%give_back_pw(v_sccs_rspace)
1273 160 : DEALLOCATE (v_sccs_rspace)
1274 : END IF
1275 :
1276 124119 : IF (gapw) THEN
1277 22362 : IF (dft_control%apply_external_potential) THEN
1278 : ! Integrals of the Hartree potential with g0_soft
1279 : CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, &
1280 42 : v_qmmm=vee, scale=-1.0_dp)
1281 : END IF
1282 22362 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces)
1283 : ! Place Vh_1c_gg_integrals after integrate_vhg0_rspace for CNEO calculations
1284 : ! because vhg0 integral is needed to build the complete nuclear equation
1285 22362 : CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c, local_rho_set=local_rho_set)
1286 : CALL Vh_1c_gg_integrals(qs_env, energy%hartree_1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE., &
1287 22362 : core_2nd=.FALSE.)
1288 : ! CNEO quantum nuclear core energy (kinetic + Z*erfc(r)/r potential from classical nuclei)
1289 22362 : energy%core_cneo = 0.0_dp
1290 22362 : IF (ASSOCIATED(local_rho_set%rhoz_cneo_set)) THEN
1291 184 : DO iatom = 1, SIZE(local_rho_set%rhoz_cneo_set)
1292 184 : energy%core_cneo = energy%core_cneo + local_rho_set%rhoz_cneo_set(iatom)%e_core
1293 : END DO
1294 : END IF
1295 : END IF
1296 :
1297 124119 : IF (gapw .OR. gapw_xc) THEN
1298 : ! Single atom contributions in the KS matrix ***
1299 26698 : CALL update_ks_atom(qs_env, ks_matrix, rho_ao, calculate_forces)
1300 26698 : IF (dft_control%do_admm) THEN
1301 : !Single atom contribution to the AUX matrices
1302 : !Note: also update ks_aux_fit matrix in case of rtp
1303 4674 : CALL admm_update_ks_atom(qs_env, calculate_forces)
1304 : END IF
1305 : END IF
1306 :
1307 : !Calculation of Mulliken restraint, if requested
1308 : CALL qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, &
1309 124119 : ks_matrix, matrix_s, rho, mulliken_order_p)
1310 :
1311 : ! Add DFT+U contribution, if requested
1312 124119 : IF (dft_control%dft_plus_u) THEN
1313 2604 : IF (just_energy) THEN
1314 746 : CALL plus_u(qs_env=qs_env)
1315 : ELSE
1316 1858 : CALL plus_u(qs_env=qs_env, matrix_h=ks_matrix)
1317 : END IF
1318 : ELSE
1319 121515 : energy%dft_plus_u = 0.0_dp
1320 : END IF
1321 :
1322 : ! At this point the ks matrix should be up to date, filter it if requested
1323 273257 : DO ispin = 1, nspins
1324 843361 : DO img = 1, nimages
1325 : CALL dbcsr_filter(ks_matrix(ispin, img)%matrix, &
1326 719242 : dft_control%qs_control%eps_filter_matrix)
1327 : END DO
1328 : END DO
1329 :
1330 : !** merge the auxiliary KS matrix and the primary one
1331 124119 : IF (dft_control%do_admm_mo) THEN
1332 13076 : IF (qs_env%run_rtp) THEN
1333 92 : CALL rtp_admm_merge_ks_matrix(qs_env)
1334 : ELSE
1335 12984 : CALL admm_mo_merge_ks_matrix(qs_env)
1336 : END IF
1337 111043 : ELSE IF (dft_control%do_admm_dm) THEN
1338 214 : CALL admm_dm_merge_ks_matrix(qs_env)
1339 : END IF
1340 :
1341 : ! External field (nonperiodic case)
1342 124119 : CALL qs_efield_local_operator(qs_env, just_energy, calculate_forces)
1343 :
1344 : ! Right now we can compute the orbital derivative here, as it depends currently only on the available
1345 : ! Kohn-Sham matrix. This might change in the future, in which case more pieces might need to be assembled
1346 : ! from this routine, notice that this part of the calculation in not linear scaling
1347 : ! right now this operation is only non-trivial because of occupation numbers and the restricted keyword
1348 124119 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy .AND. .NOT. qs_env%run_rtp) THEN
1349 45639 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs)
1350 45639 : CPASSERT(nimages == 1)
1351 45639 : ksmat => ks_matrix(:, 1)
1352 45639 : CALL calc_mo_derivatives(qs_env, ksmat, mo_derivs)
1353 : END IF
1354 :
1355 : ! ADMM overlap forces
1356 124119 : IF (calculate_forces .AND. dft_control%do_admm) THEN
1357 316 : IF (dokp) THEN
1358 30 : CALL calc_admm_ovlp_forces_kp(qs_env)
1359 : ELSE
1360 286 : CALL calc_admm_ovlp_forces(qs_env)
1361 : END IF
1362 : END IF
1363 :
1364 : ! deal with low spin roks
1365 : CALL low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
1366 124119 : calculate_forces, auxbas_pw_pool)
1367 :
1368 : ! deal with sic on explicit orbitals
1369 : CALL sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, &
1370 124119 : calculate_forces, auxbas_pw_pool)
1371 :
1372 : ! Periodic external field
1373 124119 : CALL qs_efield_berry_phase(qs_env, just_energy, calculate_forces)
1374 :
1375 : ! adds s2_restraint energy and orbital derivatives
1376 : CALL qs_ks_s2_restraint(dft_control, qs_env, matrix_s, &
1377 124119 : energy, calculate_forces, just_energy)
1378 :
1379 124119 : IF (do_ppl) THEN
1380 : ! update core energy for grid based local pseudopotential
1381 60 : ecore_ppl = 0._dp
1382 126 : DO ispin = 1, nspins
1383 126 : ecore_ppl = ecore_ppl + pw_integral_ab(vppl_rspace, rho_r(ispin))
1384 : END DO
1385 60 : energy%core = energy%core + ecore_ppl
1386 : END IF
1387 :
1388 124119 : IF (lrigpw) THEN
1389 : ! update core energy for ppl_ri method
1390 466 : CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density)
1391 466 : IF (lri_env%ppl_ri) THEN
1392 8 : ecore_ppl = 0._dp
1393 16 : DO ispin = 1, nspins
1394 8 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1395 16 : CALL v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl)
1396 : END DO
1397 8 : energy%core = energy%core + ecore_ppl
1398 : END IF
1399 : END IF
1400 :
1401 : ! Sum all energy terms to obtain the total energy
1402 : energy%total = energy%core_overlap + energy%core_self + energy%core_cneo + energy%core + &
1403 : energy%hartree + energy%hartree_1c + energy%exc + energy%exc1 + energy%ex + &
1404 : energy%dispersion + energy%gcp + energy%qmmm_el + energy%mulliken + &
1405 : SUM(energy%ddapc_restraint) + energy%s2_restraint + &
1406 : energy%dft_plus_u + energy%kTS + &
1407 : energy%efield + energy%efield_core + energy%ee + &
1408 : energy%ee_core + energy%exc_aux_fit + energy%image_charge + &
1409 248342 : energy%sccs_pol + energy%cdft + energy%exc1_aux_fit
1410 :
1411 124119 : IF (dft_control%apply_embed_pot) energy%total = energy%total + energy%embed_corr
1412 :
1413 124119 : IF (native_skala_restore_exc) energy%total = native_skala_total_scf
1414 :
1415 124119 : IF (abnormal_value(energy%total)) THEN
1416 0 : CPABORT("KS energy is an abnormal value (NaN/Inf).")
1417 : END IF
1418 :
1419 : ! Print detailed energy
1420 124119 : IF (my_print) THEN
1421 124097 : CALL print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
1422 : END IF
1423 :
1424 124119 : CALL timestop(handle)
1425 :
1426 372357 : END SUBROUTINE qs_ks_build_kohn_sham_matrix
1427 :
1428 : ! **************************************************************************************************
1429 : !> \brief ...
1430 : !> \param rho_tot_gspace ...
1431 : !> \param qs_env ...
1432 : !> \param rho ...
1433 : !> \param skip_nuclear_density ...
1434 : ! **************************************************************************************************
1435 127847 : SUBROUTINE calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
1436 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_tot_gspace
1437 : TYPE(qs_environment_type), POINTER :: qs_env
1438 : TYPE(qs_rho_type), POINTER :: rho
1439 : LOGICAL, INTENT(IN), OPTIONAL :: skip_nuclear_density
1440 :
1441 : INTEGER :: ispin
1442 : LOGICAL :: my_skip
1443 : TYPE(dft_control_type), POINTER :: dft_control
1444 127847 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
1445 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
1446 : TYPE(qs_charges_type), POINTER :: qs_charges
1447 :
1448 127847 : my_skip = .FALSE.
1449 930 : IF (PRESENT(skip_nuclear_density)) my_skip = skip_nuclear_density
1450 :
1451 127847 : CALL qs_rho_get(rho, rho_g=rho_g)
1452 127847 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
1453 :
1454 127847 : IF (.NOT. my_skip) THEN
1455 126927 : NULLIFY (rho_core)
1456 126927 : CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
1457 126927 : IF (dft_control%qs_control%gapw) THEN
1458 22682 : NULLIFY (rho0_s_gs, rhoz_cneo_s_gs)
1459 22682 : CALL get_qs_env(qs_env=qs_env, rho0_s_gs=rho0_s_gs, rhoz_cneo_s_gs=rhoz_cneo_s_gs)
1460 22682 : CPASSERT(ASSOCIATED(rho0_s_gs))
1461 22682 : CALL pw_copy(rho0_s_gs, rho_tot_gspace)
1462 22682 : IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
1463 48 : CALL pw_axpy(rhoz_cneo_s_gs, rho_tot_gspace)
1464 : END IF
1465 22682 : IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
1466 2692 : CALL pw_axpy(rho_core, rho_tot_gspace)
1467 : END IF
1468 : ELSE
1469 104245 : CALL pw_copy(rho_core, rho_tot_gspace)
1470 : END IF
1471 279183 : DO ispin = 1, dft_control%nspins
1472 279183 : CALL pw_axpy(rho_g(ispin), rho_tot_gspace)
1473 : END DO
1474 126927 : CALL get_qs_env(qs_env=qs_env, qs_charges=qs_charges)
1475 126927 : qs_charges%total_rho_gspace = pw_integrate_function(rho_tot_gspace, isign=-1)
1476 : ELSE
1477 1844 : DO ispin = 1, dft_control%nspins
1478 1844 : CALL pw_axpy(rho_g(ispin), rho_tot_gspace)
1479 : END DO
1480 : END IF
1481 :
1482 127847 : END SUBROUTINE calc_rho_tot_gspace
1483 :
1484 : ! **************************************************************************************************
1485 : !> \brief compute MO derivatives
1486 : !> \param qs_env the qs_env to update
1487 : !> \param ks_matrix ...
1488 : !> \param mo_derivs ...
1489 : !> \par History
1490 : !> 01.2014 created, transferred from qs_ks_build_kohn_sham_matrix in
1491 : !> separate subroutine
1492 : !> \author Dorothea Golze
1493 : ! **************************************************************************************************
1494 45639 : SUBROUTINE calc_mo_derivatives(qs_env, ks_matrix, mo_derivs)
1495 : TYPE(qs_environment_type), POINTER :: qs_env
1496 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, mo_derivs
1497 :
1498 : INTEGER :: ispin
1499 : LOGICAL :: uniform_occupation
1500 45639 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
1501 : TYPE(cp_fm_type), POINTER :: mo_coeff
1502 : TYPE(dbcsr_type) :: mo_derivs2_tmp1, mo_derivs2_tmp2
1503 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
1504 : TYPE(dft_control_type), POINTER :: dft_control
1505 45639 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
1506 :
1507 45639 : NULLIFY (dft_control, mo_array, mo_coeff, mo_coeff_b, occupation_numbers)
1508 :
1509 : CALL get_qs_env(qs_env, &
1510 : dft_control=dft_control, &
1511 45639 : mos=mo_array)
1512 :
1513 100071 : DO ispin = 1, SIZE(mo_derivs)
1514 :
1515 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, &
1516 54432 : mo_coeff_b=mo_coeff_b, occupation_numbers=occupation_numbers)
1517 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin)%matrix, mo_coeff_b, &
1518 54432 : 0.0_dp, mo_derivs(ispin)%matrix)
1519 :
1520 100071 : IF (dft_control%restricted) THEN
1521 : ! only the first mo_set are actual variables, but we still need both
1522 636 : CPASSERT(ispin == 1)
1523 636 : CPASSERT(SIZE(mo_array) == 2)
1524 : ! use a temporary array with the same size as the first spin for the second spin
1525 :
1526 : ! uniform_occupation is needed for this case, otherwise we can not
1527 : ! reconstruct things in ot, since we irreversibly sum
1528 636 : CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
1529 636 : CPASSERT(uniform_occupation)
1530 636 : CALL get_mo_set(mo_set=mo_array(2), uniform_occupation=uniform_occupation)
1531 636 : CPASSERT(uniform_occupation)
1532 :
1533 : ! The beta-spin might have fewer orbitals than alpa-spin...
1534 : ! create temporary matrices with beta_nmo columns
1535 636 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff_b)
1536 636 : CALL dbcsr_create(mo_derivs2_tmp1, template=mo_coeff_b)
1537 :
1538 : ! calculate beta derivatives
1539 636 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(2)%matrix, mo_coeff_b, 0.0_dp, mo_derivs2_tmp1)
1540 :
1541 : ! create larger matrix with alpha_nmo columns
1542 636 : CALL dbcsr_create(mo_derivs2_tmp2, template=mo_derivs(1)%matrix)
1543 636 : CALL dbcsr_set(mo_derivs2_tmp2, 0.0_dp)
1544 :
1545 : ! copy into larger matrix, fills the first beta_nmo columns
1546 : CALL dbcsr_copy_columns_hack(mo_derivs2_tmp2, mo_derivs2_tmp1, &
1547 : mo_array(2)%nmo, 1, 1, &
1548 : para_env=mo_array(1)%mo_coeff%matrix_struct%para_env, &
1549 636 : blacs_env=mo_array(1)%mo_coeff%matrix_struct%context)
1550 :
1551 : ! add beta contribution to alpa mo_derivs
1552 636 : CALL dbcsr_add(mo_derivs(1)%matrix, mo_derivs2_tmp2, 1.0_dp, 1.0_dp)
1553 636 : CALL dbcsr_release(mo_derivs2_tmp1)
1554 636 : CALL dbcsr_release(mo_derivs2_tmp2)
1555 : END IF
1556 : END DO
1557 :
1558 45639 : IF (dft_control%do_admm_mo) THEN
1559 6802 : CALL calc_admm_mo_derivatives(qs_env, mo_derivs)
1560 : END IF
1561 :
1562 45639 : END SUBROUTINE calc_mo_derivatives
1563 :
1564 : ! **************************************************************************************************
1565 : !> \brief updates the Kohn Sham matrix of the given qs_env (facility method)
1566 : !> \param qs_env the qs_env to update
1567 : !> \param calculate_forces if true calculate the quantities needed
1568 : !> to calculate the forces. Defaults to false.
1569 : !> \param just_energy if true updates the energies but not the
1570 : !> ks matrix. Defaults to false
1571 : !> \param print_active ...
1572 : !> \par History
1573 : !> 4.2002 created [fawzi]
1574 : !> 8.2014 kpoints [JGH]
1575 : !> 10.2014 refractored [Ole Schuett]
1576 : !> \author Fawzi Mohamed
1577 : ! **************************************************************************************************
1578 270974 : SUBROUTINE qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, &
1579 : print_active)
1580 : TYPE(qs_environment_type), POINTER :: qs_env
1581 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces, just_energy, &
1582 : print_active
1583 :
1584 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_update_qs_env'
1585 :
1586 : INTEGER :: handle, unit_nr
1587 : LOGICAL :: c_forces, do_rebuild, energy_only, &
1588 : forces_up_to_date, potential_changed, &
1589 : rho_changed, s_mstruct_changed
1590 : TYPE(qs_ks_env_type), POINTER :: ks_env
1591 :
1592 270974 : NULLIFY (ks_env)
1593 270974 : unit_nr = cp_logger_get_default_io_unit()
1594 :
1595 270974 : c_forces = .FALSE.
1596 270974 : energy_only = .FALSE.
1597 270974 : IF (PRESENT(just_energy)) energy_only = just_energy
1598 270974 : IF (PRESENT(calculate_forces)) c_forces = calculate_forces
1599 :
1600 270974 : IF (c_forces) THEN
1601 10603 : CALL timeset(routineN//'_forces', handle)
1602 : ELSE
1603 260371 : CALL timeset(routineN, handle)
1604 : END IF
1605 :
1606 270974 : CPASSERT(ASSOCIATED(qs_env))
1607 :
1608 : CALL get_qs_env(qs_env, &
1609 : ks_env=ks_env, &
1610 : rho_changed=rho_changed, &
1611 : s_mstruct_changed=s_mstruct_changed, &
1612 : potential_changed=potential_changed, &
1613 270974 : forces_up_to_date=forces_up_to_date)
1614 :
1615 270974 : do_rebuild = .FALSE.
1616 270974 : do_rebuild = do_rebuild .OR. rho_changed
1617 8801 : do_rebuild = do_rebuild .OR. s_mstruct_changed
1618 8801 : do_rebuild = do_rebuild .OR. potential_changed
1619 8801 : do_rebuild = do_rebuild .OR. (c_forces .AND. .NOT. forces_up_to_date)
1620 :
1621 : IF (do_rebuild) THEN
1622 262555 : CALL evaluate_core_matrix_traces(qs_env)
1623 :
1624 : ! the ks matrix will be rebuilt so this is fine now
1625 262555 : CALL set_ks_env(ks_env, potential_changed=.FALSE.)
1626 :
1627 : CALL rebuild_ks_matrix(qs_env, &
1628 : calculate_forces=c_forces, &
1629 : just_energy=energy_only, &
1630 262555 : print_active=print_active)
1631 :
1632 262555 : IF (.NOT. energy_only) THEN
1633 : CALL set_ks_env(ks_env, &
1634 : rho_changed=.FALSE., &
1635 : s_mstruct_changed=.FALSE., &
1636 484455 : forces_up_to_date=forces_up_to_date .OR. c_forces)
1637 : END IF
1638 : END IF
1639 :
1640 270974 : CALL timestop(handle)
1641 :
1642 270974 : END SUBROUTINE qs_ks_update_qs_env
1643 :
1644 : ! **************************************************************************************************
1645 : !> \brief Calculates the traces of the core matrices and the density matrix.
1646 : !> \param qs_env ...
1647 : !> \param rho_ao_ext ...
1648 : !> \author Ole Schuett
1649 : ! **************************************************************************************************
1650 289687 : SUBROUTINE evaluate_core_matrix_traces(qs_env, rho_ao_ext)
1651 : TYPE(qs_environment_type), POINTER :: qs_env
1652 : TYPE(dbcsr_p_type), DIMENSION(:, :), OPTIONAL, &
1653 : POINTER :: rho_ao_ext
1654 :
1655 : CHARACTER(LEN=*), PARAMETER :: routineN = 'evaluate_core_matrix_traces'
1656 :
1657 : INTEGER :: handle
1658 : REAL(KIND=dp) :: energy_core_im
1659 289687 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_t, rho_ao_kp
1660 : TYPE(dft_control_type), POINTER :: dft_control
1661 : TYPE(qs_energy_type), POINTER :: energy
1662 : TYPE(qs_rho_type), POINTER :: rho
1663 :
1664 289687 : CALL timeset(routineN, handle)
1665 289687 : NULLIFY (energy, rho, dft_control, rho_ao_kp, matrixkp_t, matrixkp_h)
1666 :
1667 : CALL get_qs_env(qs_env, &
1668 : rho=rho, &
1669 : energy=energy, &
1670 : dft_control=dft_control, &
1671 : kinetic_kp=matrixkp_t, &
1672 289687 : matrix_h_kp=matrixkp_h)
1673 :
1674 289687 : IF (PRESENT(rho_ao_ext)) THEN
1675 27056 : rho_ao_kp => rho_ao_ext
1676 : ELSE
1677 262631 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1678 : END IF
1679 :
1680 289687 : CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy%core, dft_control%nspins)
1681 :
1682 : ! Add the imaginary part in the RTP case
1683 289687 : IF (qs_env%run_rtp) THEN
1684 3220 : IF (dft_control%rtp_control%velocity_gauge) THEN
1685 150 : CALL get_qs_env(qs_env, matrix_h_im_kp=matrixkp_h)
1686 150 : CALL qs_rho_get(rho, rho_ao_im_kp=rho_ao_kp)
1687 150 : CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy_core_im, dft_control%nspins)
1688 150 : energy%core = energy%core - energy_core_im
1689 : END IF
1690 : END IF
1691 :
1692 : ! kinetic energy
1693 289687 : IF (ASSOCIATED(matrixkp_t)) THEN
1694 123809 : CALL calculate_ptrace(matrixkp_t, rho_ao_kp, energy%kinetic, dft_control%nspins)
1695 : END IF
1696 :
1697 289687 : CALL timestop(handle)
1698 289687 : END SUBROUTINE evaluate_core_matrix_traces
1699 :
1700 : ! **************************************************************************************************
1701 : !> \brief Constructs a new Khon-Sham matrix
1702 : !> \param qs_env ...
1703 : !> \param calculate_forces ...
1704 : !> \param just_energy ...
1705 : !> \param print_active ...
1706 : !> \author Ole Schuett
1707 : ! **************************************************************************************************
1708 262575 : SUBROUTINE rebuild_ks_matrix(qs_env, calculate_forces, just_energy, print_active)
1709 : TYPE(qs_environment_type), POINTER :: qs_env
1710 : LOGICAL, INTENT(IN) :: calculate_forces, just_energy
1711 : LOGICAL, INTENT(IN), OPTIONAL :: print_active
1712 :
1713 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rebuild_ks_matrix'
1714 :
1715 : INTEGER :: handle
1716 : TYPE(dft_control_type), POINTER :: dft_control
1717 :
1718 262575 : CALL timeset(routineN, handle)
1719 262575 : NULLIFY (dft_control)
1720 :
1721 262575 : CALL get_qs_env(qs_env, dft_control=dft_control)
1722 :
1723 262575 : IF (dft_control%qs_control%semi_empirical) THEN
1724 : CALL build_se_fock_matrix(qs_env, &
1725 : calculate_forces=calculate_forces, &
1726 41348 : just_energy=just_energy)
1727 :
1728 221227 : ELSE IF (dft_control%qs_control%dftb) THEN
1729 : CALL build_dftb_ks_matrix(qs_env, &
1730 : calculate_forces=calculate_forces, &
1731 30416 : just_energy=just_energy)
1732 :
1733 190811 : ELSE IF (dft_control%qs_control%xtb) THEN
1734 66990 : IF (dft_control%qs_control%xtb_control%do_tblite) THEN
1735 : CALL build_tblite_ks_matrix(qs_env, &
1736 : calculate_forces=calculate_forces, &
1737 29808 : just_energy=just_energy)
1738 : ELSE
1739 : CALL build_xtb_ks_matrix(qs_env, &
1740 : calculate_forces=calculate_forces, &
1741 37182 : just_energy=just_energy)
1742 : END IF
1743 : ELSE
1744 : CALL qs_ks_build_kohn_sham_matrix(qs_env, &
1745 : calculate_forces=calculate_forces, &
1746 : just_energy=just_energy, &
1747 123821 : print_active=print_active)
1748 : END IF
1749 :
1750 262575 : CALL timestop(handle)
1751 :
1752 262575 : END SUBROUTINE rebuild_ks_matrix
1753 :
1754 : ! **************************************************************************************************
1755 : !> \brief Allocate ks_matrix if necessary, take current overlap matrix as template
1756 : !> \param qs_env ...
1757 : !> \param is_complex ...
1758 : !> \par History
1759 : !> refactoring 04.03.2011 [MI]
1760 : !> \author
1761 : ! **************************************************************************************************
1762 :
1763 28706 : SUBROUTINE qs_ks_allocate_basics(qs_env, is_complex)
1764 : TYPE(qs_environment_type), POINTER :: qs_env
1765 : LOGICAL, INTENT(in) :: is_complex
1766 :
1767 : CHARACTER(LEN=default_string_length) :: headline
1768 : INTEGER :: ic, ispin, nimages, nspins
1769 : LOGICAL :: do_kpoints
1770 28706 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrixkp_im_ks, matrixkp_ks
1771 : TYPE(dbcsr_type), POINTER :: refmatrix
1772 : TYPE(dft_control_type), POINTER :: dft_control
1773 : TYPE(kpoint_type), POINTER :: kpoints
1774 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1775 28706 : POINTER :: sab_orb
1776 : TYPE(qs_ks_env_type), POINTER :: ks_env
1777 :
1778 28706 : NULLIFY (dft_control, ks_env, matrix_s_kp, sab_orb, matrixkp_ks, refmatrix, matrixkp_im_ks, kpoints)
1779 :
1780 : CALL get_qs_env(qs_env, &
1781 : dft_control=dft_control, &
1782 : matrix_s_kp=matrix_s_kp, &
1783 : ks_env=ks_env, &
1784 : kpoints=kpoints, &
1785 : do_kpoints=do_kpoints, &
1786 : matrix_ks_kp=matrixkp_ks, &
1787 28706 : matrix_ks_im_kp=matrixkp_im_ks)
1788 :
1789 28706 : IF (do_kpoints) THEN
1790 3562 : CALL get_kpoint_info(kpoints, sab_nl=sab_orb)
1791 : ELSE
1792 25144 : CALL get_qs_env(qs_env, sab_orb=sab_orb)
1793 : END IF
1794 :
1795 28706 : nspins = dft_control%nspins
1796 28706 : nimages = dft_control%nimages
1797 :
1798 28706 : IF (.NOT. ASSOCIATED(matrixkp_ks)) THEN
1799 28664 : CALL dbcsr_allocate_matrix_set(matrixkp_ks, nspins, nimages)
1800 28664 : refmatrix => matrix_s_kp(1, 1)%matrix
1801 61170 : DO ispin = 1, nspins
1802 314188 : DO ic = 1, nimages
1803 253018 : IF (nspins > 1) THEN
1804 38516 : IF (ispin == 1) THEN
1805 19258 : headline = "KOHN-SHAM MATRIX FOR ALPHA SPIN"
1806 : ELSE
1807 19258 : headline = "KOHN-SHAM MATRIX FOR BETA SPIN"
1808 : END IF
1809 : ELSE
1810 214502 : headline = "KOHN-SHAM MATRIX"
1811 : END IF
1812 253018 : ALLOCATE (matrixkp_ks(ispin, ic)%matrix)
1813 : CALL dbcsr_create(matrix=matrixkp_ks(ispin, ic)%matrix, template=refmatrix, &
1814 253018 : name=TRIM(headline), matrix_type=dbcsr_type_symmetric)
1815 253018 : CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_ks(ispin, ic)%matrix, sab_orb)
1816 285524 : CALL dbcsr_set(matrixkp_ks(ispin, ic)%matrix, 0.0_dp)
1817 : END DO
1818 : END DO
1819 28664 : CALL set_ks_env(ks_env, matrix_ks_kp=matrixkp_ks)
1820 : END IF
1821 :
1822 28706 : IF (is_complex) THEN
1823 144 : IF (.NOT. ASSOCIATED(matrixkp_im_ks)) THEN
1824 144 : CPASSERT(nspins == SIZE(matrixkp_ks, 1))
1825 144 : CPASSERT(nimages == SIZE(matrixkp_ks, 2))
1826 144 : CALL dbcsr_allocate_matrix_set(matrixkp_im_ks, nspins, nimages)
1827 306 : DO ispin = 1, nspins
1828 468 : DO ic = 1, nimages
1829 162 : IF (nspins > 1) THEN
1830 36 : IF (ispin == 1) THEN
1831 18 : headline = "IMAGINARY KOHN-SHAM MATRIX FOR ALPHA SPIN"
1832 : ELSE
1833 18 : headline = "IMAGINARY KOHN-SHAM MATRIX FOR BETA SPIN"
1834 : END IF
1835 : ELSE
1836 126 : headline = "IMAGINARY KOHN-SHAM MATRIX"
1837 : END IF
1838 162 : ALLOCATE (matrixkp_im_ks(ispin, ic)%matrix)
1839 162 : refmatrix => matrixkp_ks(ispin, ic)%matrix ! base on real part, but anti-symmetric
1840 : CALL dbcsr_create(matrix=matrixkp_im_ks(ispin, ic)%matrix, template=refmatrix, &
1841 162 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric)
1842 162 : CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_im_ks(ispin, ic)%matrix, sab_orb)
1843 324 : CALL dbcsr_set(matrixkp_im_ks(ispin, ic)%matrix, 0.0_dp)
1844 : END DO
1845 : END DO
1846 144 : CALL set_ks_env(ks_env, matrix_ks_im_kp=matrixkp_im_ks)
1847 : END IF
1848 : END IF
1849 :
1850 28706 : END SUBROUTINE qs_ks_allocate_basics
1851 :
1852 : END MODULE qs_ks_methods
|
