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