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