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