Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 admm_dm_methods, ONLY: admm_dm_calc_rho_aux,&
24 : admm_dm_merge_ks_matrix
25 : USE admm_methods, ONLY: admm_mo_calc_rho_aux,&
26 : admm_mo_calc_rho_aux_kp,&
27 : admm_mo_merge_ks_matrix,&
28 : admm_update_ks_atom,&
29 : calc_admm_mo_derivatives,&
30 : calc_admm_ovlp_forces,&
31 : calc_admm_ovlp_forces_kp
32 : USE admm_types, ONLY: admm_type,&
33 : get_admm_env
34 : USE cell_types, ONLY: cell_type
35 : USE cp_control_types, ONLY: dft_control_type
36 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
37 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
38 : dbcsr_copy_columns_hack
39 : USE cp_ddapc, ONLY: qs_ks_ddapc
40 : USE cp_fm_types, ONLY: cp_fm_type
41 : USE cp_log_handling, ONLY: cp_get_default_logger,&
42 : cp_logger_get_default_io_unit,&
43 : cp_logger_type
44 : USE cp_output_handling, ONLY: cp_p_file,&
45 : cp_print_key_should_output
46 : USE dbcsr_api, ONLY: &
47 : dbcsr_add, dbcsr_copy, dbcsr_create, dbcsr_filter, dbcsr_get_info, dbcsr_multiply, &
48 : dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_antisymmetric, &
49 : dbcsr_type_symmetric
50 : USE dft_plus_u, ONLY: plus_u
51 : USE hartree_local_methods, ONLY: Vh_1c_gg_integrals
52 : USE hartree_local_types, ONLY: ecoul_1center_type
53 : USE hfx_admm_utils, ONLY: hfx_admm_init,&
54 : hfx_ks_matrix,&
55 : hfx_ks_matrix_kp
56 : USE input_constants, ONLY: do_ppl_grid,&
57 : outer_scf_becke_constraint,&
58 : outer_scf_hirshfeld_constraint
59 : USE input_section_types, ONLY: section_vals_get,&
60 : section_vals_get_subs_vals,&
61 : section_vals_type,&
62 : section_vals_val_get
63 : USE kg_correction, ONLY: kg_ekin_subset
64 : USE kinds, ONLY: default_string_length,&
65 : dp
66 : USE kpoint_types, ONLY: get_kpoint_info,&
67 : kpoint_type
68 : USE lri_environment_methods, ONLY: v_int_ppl_energy
69 : USE lri_environment_types, ONLY: lri_density_type,&
70 : lri_environment_type,&
71 : lri_kind_type
72 : USE mathlib, ONLY: abnormal_value
73 : USE message_passing, ONLY: mp_para_env_type
74 : USE pw_env_types, ONLY: pw_env_get,&
75 : pw_env_type
76 : USE pw_methods, ONLY: pw_axpy,&
77 : pw_copy,&
78 : pw_integral_ab,&
79 : pw_integrate_function,&
80 : pw_scale,&
81 : pw_transfer,&
82 : pw_zero
83 : USE pw_poisson_methods, ONLY: pw_poisson_solve
84 : USE pw_poisson_types, ONLY: pw_poisson_implicit,&
85 : pw_poisson_type
86 : USE pw_pool_types, ONLY: pw_pool_type
87 : USE pw_types, ONLY: pw_c1d_gs_type,&
88 : pw_r3d_rs_type
89 : USE qmmm_image_charge, ONLY: add_image_pot_to_hartree_pot,&
90 : calculate_image_pot,&
91 : integrate_potential_devga_rspace
92 : USE qs_cdft_types, ONLY: cdft_control_type
93 : USE qs_charges_types, ONLY: qs_charges_type
94 : USE qs_core_energies, ONLY: calculate_ptrace
95 : USE qs_dftb_matrices, ONLY: build_dftb_ks_matrix
96 : USE qs_efield_berry, ONLY: qs_efield_berry_phase
97 : USE qs_efield_local, ONLY: qs_efield_local_operator
98 : USE qs_energy_types, ONLY: qs_energy_type
99 : USE qs_environment_types, ONLY: get_qs_env,&
100 : qs_environment_type
101 : USE qs_gapw_densities, ONLY: prepare_gapw_den
102 : USE qs_integrate_potential, ONLY: integrate_ppl_rspace,&
103 : integrate_rho_nlcc,&
104 : integrate_v_core_rspace
105 : USE qs_ks_apply_restraints, ONLY: qs_ks_cdft_constraint,&
106 : qs_ks_mulliken_restraint,&
107 : qs_ks_s2_restraint
108 : USE qs_ks_atom, ONLY: update_ks_atom
109 : USE qs_ks_qmmm_methods, ONLY: qmmm_calculate_energy,&
110 : qmmm_modify_hartree_pot
111 : USE qs_ks_types, ONLY: qs_ks_env_type,&
112 : set_ks_env
113 : USE qs_ks_utils, ONLY: &
114 : calc_v_sic_rspace, calculate_zmp_potential, compute_matrix_vxc, &
115 : get_embed_potential_energy, low_spin_roks, print_densities, print_detailed_energy, &
116 : sic_explicit_orbitals, sum_up_and_integrate
117 : USE qs_local_rho_types, ONLY: local_rho_type
118 : USE qs_mo_types, ONLY: get_mo_set,&
119 : mo_set_type
120 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
121 : USE qs_rho0_ggrid, ONLY: integrate_vhg0_rspace
122 : USE qs_rho_types, ONLY: qs_rho_get,&
123 : qs_rho_type
124 : USE qs_sccs, ONLY: sccs
125 : USE qs_vxc, ONLY: qs_vxc_create
126 : USE qs_vxc_atom, ONLY: calculate_vxc_atom
127 : USE rtp_admm_methods, ONLY: rtp_admm_calc_rho_aux,&
128 : rtp_admm_merge_ks_matrix
129 : USE se_fock_matrix, ONLY: build_se_fock_matrix
130 : USE surface_dipole, ONLY: calc_dipsurf_potential
131 : USE virial_types, ONLY: virial_type
132 : USE xtb_matrices, ONLY: build_xtb_ks_matrix
133 : #include "./base/base_uses.f90"
134 :
135 : IMPLICIT NONE
136 :
137 : PRIVATE
138 :
139 : LOGICAL, PARAMETER :: debug_this_module = .TRUE.
140 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_methods'
141 :
142 : PUBLIC :: calc_rho_tot_gspace, qs_ks_update_qs_env, qs_ks_build_kohn_sham_matrix, &
143 : qs_ks_allocate_basics
144 :
145 : CONTAINS
146 :
147 : ! **************************************************************************************************
148 : !> \brief routine where the real calculations are made: the
149 : !> KS matrix is calculated
150 : !> \param qs_env the qs_env to update
151 : !> \param calculate_forces if true calculate the quantities needed
152 : !> to calculate the forces. Defaults to false.
153 : !> \param just_energy if true updates the energies but not the
154 : !> ks matrix. Defaults to false
155 : !> \param print_active ...
156 : !> \param ext_ks_matrix ...
157 : !> \par History
158 : !> 06.2002 moved from qs_scf to qs_ks_methods, use of ks_env
159 : !> new did_change scheme [fawzi]
160 : !> 10.2002 introduced pools, uses updated rho as input, LSD [fawzi]
161 : !> 10.2004 build_kohn_sham matrix now also computes the derivatives
162 : !> of the total energy wrt to the MO coefs, if instructed to
163 : !> do so. This appears useful for orbital dependent functionals
164 : !> where the KS matrix alone (however this might be defined)
165 : !> does not contain the info to construct this derivative.
166 : !> \author Matthias Krack
167 : !> \note
168 : !> make rho, energy and qs_charges optional, defaulting
169 : !> to qs_env components?
170 : ! **************************************************************************************************
171 289353 : SUBROUTINE qs_ks_build_kohn_sham_matrix(qs_env, calculate_forces, just_energy, &
172 : print_active, ext_ks_matrix)
173 : TYPE(qs_environment_type), POINTER :: qs_env
174 : LOGICAL, INTENT(in) :: calculate_forces, just_energy
175 : LOGICAL, INTENT(IN), OPTIONAL :: print_active
176 : TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
177 : POINTER :: ext_ks_matrix
178 :
179 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_build_kohn_sham_matrix'
180 :
181 : CHARACTER(len=default_string_length) :: name
182 : INTEGER :: handle, iatom, img, ispin, nimages, ns, &
183 : nspins
184 : LOGICAL :: do_adiabatic_rescaling, do_ddapc, do_hfx, do_ppl, dokp, gapw, gapw_xc, &
185 : hfx_treat_lsd_in_core, just_energy_xc, lrigpw, my_print, rigpw, use_virial
186 : REAL(KIND=dp) :: ecore_ppl, edisp, ee_ener, ekin_mol, &
187 : mulliken_order_p, vscale
188 : REAL(KIND=dp), DIMENSION(3, 3) :: h_stress, pv_loc
189 : TYPE(admm_type), POINTER :: admm_env
190 : TYPE(cdft_control_type), POINTER :: cdft_control
191 : TYPE(cell_type), POINTER :: cell
192 : TYPE(cp_logger_type), POINTER :: logger
193 96451 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, matrix_vxc, mo_derivs
194 96451 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix, ks_matrix_im, matrix_h, &
195 96451 : matrix_h_im, matrix_s, my_rho, rho_ao
196 : TYPE(dft_control_type), POINTER :: dft_control
197 96451 : TYPE(ecoul_1center_type), DIMENSION(:), POINTER :: ecoul_1c
198 : TYPE(local_rho_type), POINTER :: local_rho_set
199 : TYPE(lri_density_type), POINTER :: lri_density
200 : TYPE(lri_environment_type), POINTER :: lri_env
201 96451 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
202 : TYPE(mp_para_env_type), POINTER :: para_env
203 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_hartree_gspace
204 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
205 : TYPE(pw_env_type), POINTER :: pw_env
206 : TYPE(pw_poisson_type), POINTER :: poisson_env
207 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
208 96451 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace_embed, v_rspace_new, &
209 96451 : v_rspace_new_aux_fit, v_tau_rspace, &
210 96451 : v_tau_rspace_aux_fit
211 : TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs, rho_nlcc, v_hartree_rspace, &
212 : v_sccs_rspace, v_sic_rspace, &
213 : v_spin_ddapc_rest_r, vee, vppl_rspace
214 : TYPE(qs_energy_type), POINTER :: energy
215 : TYPE(qs_ks_env_type), POINTER :: ks_env
216 : TYPE(qs_rho_type), POINTER :: rho, rho_struct, rho_xc
217 : TYPE(section_vals_type), POINTER :: adiabatic_rescaling_section, &
218 : hfx_sections, input, scf_section, &
219 : xc_section
220 : TYPE(virial_type), POINTER :: virial
221 :
222 96451 : CALL timeset(routineN, handle)
223 96451 : NULLIFY (admm_env, cell, dft_control, logger, mo_derivs, my_rho, &
224 96451 : rho_struct, para_env, pw_env, virial, vppl_rspace, &
225 96451 : adiabatic_rescaling_section, hfx_sections, &
226 96451 : input, scf_section, xc_section, matrix_h, matrix_h_im, matrix_s, &
227 96451 : auxbas_pw_pool, poisson_env, v_rspace_new, v_rspace_new_aux_fit, &
228 96451 : v_tau_rspace, v_tau_rspace_aux_fit, matrix_vxc, vee, rho_nlcc, &
229 96451 : ks_env, ks_matrix, ks_matrix_im, rho, energy, rho_xc, rho_r, rho_ao, rho_core)
230 :
231 96451 : CPASSERT(ASSOCIATED(qs_env))
232 :
233 96451 : logger => cp_get_default_logger()
234 96451 : my_print = .TRUE.
235 96451 : IF (PRESENT(print_active)) my_print = print_active
236 :
237 : CALL get_qs_env(qs_env, &
238 : ks_env=ks_env, &
239 : dft_control=dft_control, &
240 : matrix_h_kp=matrix_h, &
241 : matrix_h_im_kp=matrix_h_im, &
242 : matrix_s_kp=matrix_s, &
243 : matrix_ks_kp=ks_matrix, &
244 : matrix_ks_im_kp=ks_matrix_im, &
245 : matrix_vxc=matrix_vxc, &
246 : pw_env=pw_env, &
247 : cell=cell, &
248 : para_env=para_env, &
249 : input=input, &
250 : virial=virial, &
251 : v_hartree_rspace=v_hartree_rspace, &
252 : vee=vee, &
253 : rho_nlcc=rho_nlcc, &
254 : rho=rho, &
255 : rho_core=rho_core, &
256 : rho_xc=rho_xc, &
257 96451 : energy=energy)
258 :
259 96451 : CALL qs_rho_get(rho, rho_r=rho_r, rho_ao_kp=rho_ao)
260 :
261 96451 : IF (PRESENT(ext_ks_matrix)) THEN
262 : ! remap pointer to allow for non-kpoint external ks matrix
263 : ! ext_ks_matrix is used in linear response code
264 32 : ns = SIZE(ext_ks_matrix)
265 32 : ks_matrix(1:ns, 1:1) => ext_ks_matrix(1:ns)
266 : END IF
267 :
268 96451 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
269 :
270 96451 : hfx_sections => section_vals_get_subs_vals(input, "DFT%XC%HF")
271 96451 : CALL section_vals_get(hfx_sections, explicit=do_hfx)
272 96451 : IF (do_hfx) THEN
273 : CALL section_vals_val_get(hfx_sections, "TREAT_LSD_IN_CORE", l_val=hfx_treat_lsd_in_core, &
274 23518 : i_rep_section=1)
275 : END IF
276 96451 : adiabatic_rescaling_section => section_vals_get_subs_vals(input, "DFT%XC%ADIABATIC_RESCALING")
277 96451 : CALL section_vals_get(adiabatic_rescaling_section, explicit=do_adiabatic_rescaling)
278 96451 : just_energy_xc = just_energy
279 96451 : IF (do_adiabatic_rescaling) THEN
280 : !! If we perform adiabatic rescaling, the xc potential has to be scaled by the xc- and
281 : !! HFX-energy. Thus, let us first calculate the energy
282 56 : just_energy_xc = .TRUE.
283 : END IF
284 :
285 96451 : nimages = dft_control%nimages
286 96451 : nspins = dft_control%nspins
287 96451 : CPASSERT(ASSOCIATED(matrix_h))
288 96451 : CPASSERT(ASSOCIATED(matrix_s))
289 96451 : CPASSERT(ASSOCIATED(rho))
290 96451 : CPASSERT(ASSOCIATED(pw_env))
291 96451 : CPASSERT(SIZE(ks_matrix, 1) > 0)
292 96451 : dokp = (nimages > 1)
293 :
294 : ! Setup the possible usage of DDAPC charges
295 : do_ddapc = dft_control%qs_control%ddapc_restraint .OR. &
296 : qs_env%cp_ddapc_ewald%do_decoupling .OR. &
297 : qs_env%cp_ddapc_ewald%do_qmmm_periodic_decpl .OR. &
298 96451 : qs_env%cp_ddapc_ewald%do_solvation
299 :
300 : ! Check if LRIGPW is used
301 96451 : lrigpw = dft_control%qs_control%lrigpw
302 96451 : rigpw = dft_control%qs_control%rigpw
303 96451 : IF (rigpw) THEN
304 0 : CPASSERT(nimages == 1)
305 : END IF
306 96451 : IF (lrigpw .AND. rigpw) THEN
307 0 : CPABORT(" LRI and RI are not compatible")
308 : END IF
309 :
310 : ! Check for GAPW method : additional terms for local densities
311 96451 : gapw = dft_control%qs_control%gapw
312 96451 : gapw_xc = dft_control%qs_control%gapw_xc
313 96451 : IF (gapw_xc .AND. gapw) THEN
314 0 : CPABORT(" GAPW and GAPW_XC are not compatible")
315 : END IF
316 96451 : IF ((gapw .AND. lrigpw) .OR. (gapw_xc .AND. lrigpw)) THEN
317 0 : CPABORT(" GAPW/GAPW_XC and LRIGPW are not compatible")
318 : END IF
319 96451 : IF ((gapw .AND. rigpw) .OR. (gapw_xc .AND. rigpw)) THEN
320 0 : CPABORT(" GAPW/GAPW_XC and RIGPW are not compatible")
321 : END IF
322 :
323 96451 : do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
324 96451 : IF (do_ppl) THEN
325 60 : CPASSERT(.NOT. gapw)
326 60 : CALL get_qs_env(qs_env=qs_env, vppl=vppl_rspace)
327 : END IF
328 :
329 96451 : IF (gapw_xc) THEN
330 2534 : CPASSERT(ASSOCIATED(rho_xc))
331 : END IF
332 :
333 : ! gets the tmp grids
334 96451 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env)
335 :
336 96451 : IF (gapw .AND. (poisson_env%parameters%solver .EQ. pw_poisson_implicit)) THEN
337 0 : CPABORT("The implicit Poisson solver cannot be used in conjunction with GAPW.")
338 : END IF
339 :
340 : ! *** Prepare densities for gapw ***
341 96451 : IF (gapw .OR. gapw_xc) THEN
342 15506 : CALL prepare_gapw_den(qs_env, do_rho0=(.NOT. gapw_xc))
343 : END IF
344 :
345 : ! Calculate the Hartree potential
346 96451 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
347 96451 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
348 :
349 96451 : scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
350 : IF (BTEST(cp_print_key_should_output(logger%iter_info, scf_section, &
351 : "PRINT%DETAILED_ENERGY"), &
352 : cp_p_file) .AND. &
353 96451 : (.NOT. gapw) .AND. (.NOT. gapw_xc) .AND. &
354 : (.NOT. (poisson_env%parameters%solver .EQ. pw_poisson_implicit))) THEN
355 904 : CALL pw_zero(rho_tot_gspace)
356 904 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density=.TRUE.)
357 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%e_hartree, &
358 904 : v_hartree_gspace)
359 904 : CALL pw_zero(rho_tot_gspace)
360 904 : CALL pw_zero(v_hartree_gspace)
361 : END IF
362 :
363 : ! Get the total density in g-space [ions + electrons]
364 96451 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
365 :
366 96451 : IF (my_print) THEN
367 96429 : CALL print_densities(qs_env, rho)
368 : END IF
369 :
370 96451 : IF (dft_control%do_sccs) THEN
371 : ! Self-consistent continuum solvation (SCCS) model
372 : NULLIFY (v_sccs_rspace)
373 132 : ALLOCATE (v_sccs_rspace)
374 132 : CALL auxbas_pw_pool%create_pw(v_sccs_rspace)
375 :
376 132 : IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
377 0 : CPABORT("The implicit Poisson solver cannot be used together with SCCS.")
378 : END IF
379 :
380 132 : IF (use_virial .AND. calculate_forces) THEN
381 : CALL sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs_rspace, &
382 0 : h_stress=h_stress)
383 0 : virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
384 0 : virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
385 : ELSE
386 132 : CALL sccs(qs_env, rho_tot_gspace, v_hartree_gspace, v_sccs_rspace)
387 : END IF
388 : ELSE
389 : ! Getting the Hartree energy and Hartree potential. Also getting the stress tensor
390 : ! from the Hartree term if needed. No nuclear force information here
391 96319 : IF (use_virial .AND. calculate_forces) THEN
392 374 : h_stress(:, :) = 0.0_dp
393 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%hartree, &
394 : v_hartree_gspace, h_stress=h_stress, &
395 374 : rho_core=rho_core)
396 4862 : virial%pv_ehartree = virial%pv_ehartree + h_stress/REAL(para_env%num_pe, dp)
397 4862 : virial%pv_virial = virial%pv_virial + h_stress/REAL(para_env%num_pe, dp)
398 : ELSE
399 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, energy%hartree, &
400 95945 : v_hartree_gspace, rho_core=rho_core)
401 : END IF
402 : END IF
403 :
404 : ! In case decouple periodic images and/or apply restraints to charges
405 96451 : IF (do_ddapc) THEN
406 : CALL qs_ks_ddapc(qs_env, auxbas_pw_pool, rho_tot_gspace, v_hartree_gspace, &
407 : v_spin_ddapc_rest_r, energy, calculate_forces, ks_matrix, &
408 1016 : just_energy)
409 : ELSE
410 95435 : dft_control%qs_control%ddapc_explicit_potential = .FALSE.
411 95435 : dft_control%qs_control%ddapc_restraint_is_spin = .FALSE.
412 95435 : IF (.NOT. just_energy) THEN
413 88718 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
414 88718 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
415 : END IF
416 : END IF
417 96451 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
418 :
419 96451 : IF (dft_control%correct_surf_dip) THEN
420 98 : IF (dft_control%surf_dip_correct_switch) THEN
421 98 : CALL calc_dipsurf_potential(qs_env, energy)
422 98 : energy%hartree = energy%hartree + energy%surf_dipole
423 : END IF
424 : END IF
425 :
426 : ! SIC
427 : CALL calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, &
428 96451 : just_energy, calculate_forces, auxbas_pw_pool)
429 :
430 96451 : IF (gapw) THEN
431 12972 : CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c, local_rho_set=local_rho_set)
432 : CALL Vh_1c_gg_integrals(qs_env, energy%hartree_1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE., &
433 12972 : core_2nd=.FALSE.)
434 : END IF
435 :
436 : ! Check if CDFT constraint is needed
437 96451 : CALL qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
438 :
439 : ! Adds the External Potential if requested
440 96451 : IF (dft_control%apply_external_potential) THEN
441 : ! Compute the energy due to the external potential
442 : ee_ener = 0.0_dp
443 728 : DO ispin = 1, nspins
444 728 : ee_ener = ee_ener + pw_integral_ab(rho_r(ispin), vee)
445 : END DO
446 364 : IF (.NOT. just_energy) THEN
447 364 : IF (gapw) THEN
448 : CALL get_qs_env(qs_env=qs_env, &
449 42 : rho0_s_rs=rho0_s_rs)
450 42 : CPASSERT(ASSOCIATED(rho0_s_rs))
451 42 : ee_ener = ee_ener + pw_integral_ab(rho0_s_rs, vee)
452 : END IF
453 : END IF
454 : ! the sign accounts for the charge of the electrons
455 364 : energy%ee = -ee_ener
456 : END IF
457 :
458 : ! Adds the QM/MM potential
459 96451 : IF (qs_env%qmmm) THEN
460 : CALL qmmm_calculate_energy(qs_env=qs_env, &
461 : rho=rho_r, &
462 : v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, &
463 6298 : qmmm_energy=energy%qmmm_el)
464 6298 : IF (qs_env%qmmm_env_qm%image_charge) THEN
465 : CALL calculate_image_pot(v_hartree_rspace=v_hartree_rspace, &
466 : rho_hartree_gspace=rho_tot_gspace, &
467 : energy=energy, &
468 : qmmm_env=qs_env%qmmm_env_qm, &
469 60 : qs_env=qs_env)
470 60 : IF (.NOT. just_energy) THEN
471 : CALL add_image_pot_to_hartree_pot(v_hartree=v_hartree_rspace, &
472 : v_metal=qs_env%ks_qmmm_env%v_metal_rspace, &
473 60 : qs_env=qs_env)
474 60 : IF (calculate_forces) THEN
475 : CALL integrate_potential_devga_rspace( &
476 : potential=v_hartree_rspace, coeff=qs_env%image_coeff, &
477 : forces=qs_env%qmmm_env_qm%image_charge_pot%image_forcesMM, &
478 20 : qmmm_env=qs_env%qmmm_env_qm, qs_env=qs_env)
479 : END IF
480 : END IF
481 60 : CALL qs_env%ks_qmmm_env%v_metal_rspace%release()
482 60 : DEALLOCATE (qs_env%ks_qmmm_env%v_metal_rspace)
483 : END IF
484 6298 : IF (.NOT. just_energy) THEN
485 : CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, &
486 6218 : v_qmmm=qs_env%ks_qmmm_env%v_qmmm_rspace, scale=1.0_dp)
487 : END IF
488 : END IF
489 96451 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
490 :
491 : ! calculate the density matrix for the fitted mo_coeffs
492 96451 : IF (dft_control%do_admm) THEN
493 10180 : CALL hfx_admm_init(qs_env, calculate_forces)
494 :
495 10180 : IF (dft_control%do_admm_mo) THEN
496 10040 : IF (qs_env%run_rtp) THEN
497 76 : CALL rtp_admm_calc_rho_aux(qs_env)
498 : ELSE
499 9964 : IF (dokp) THEN
500 114 : CALL admm_mo_calc_rho_aux_kp(qs_env)
501 : ELSE
502 9850 : CALL admm_mo_calc_rho_aux(qs_env)
503 : END IF
504 : END IF
505 140 : ELSEIF (dft_control%do_admm_dm) THEN
506 140 : CALL admm_dm_calc_rho_aux(qs_env)
507 : END IF
508 : END IF
509 :
510 : ! only activate stress calculation if
511 96451 : IF (use_virial .AND. calculate_forces) virial%pv_calculate = .TRUE.
512 :
513 : ! *** calculate the xc potential on the pw density ***
514 : ! *** associates v_rspace_new if the xc potential needs to be computed.
515 : ! If we do wavefunction fitting, we need the vxc_potential in the auxiliary basis set
516 96451 : IF (dft_control%do_admm) THEN
517 10180 : CALL get_qs_env(qs_env, admm_env=admm_env)
518 10180 : xc_section => admm_env%xc_section_aux
519 10180 : CALL get_admm_env(admm_env, rho_aux_fit=rho_struct)
520 :
521 : ! here we ignore a possible vdW section in admm_env%xc_section_aux
522 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
523 : vxc_rho=v_rspace_new_aux_fit, vxc_tau=v_tau_rspace_aux_fit, exc=energy%exc_aux_fit, &
524 10180 : just_energy=just_energy_xc)
525 :
526 10180 : IF (admm_env%do_gapw) THEN
527 : !compute the potential due to atomic densities
528 : CALL calculate_vxc_atom(qs_env, energy_only=just_energy_xc, exc1=energy%exc1_aux_fit, &
529 : kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
530 : xc_section_external=xc_section, &
531 2350 : rho_atom_set_external=admm_env%admm_gapw_env%local_rho_set%rho_atom_set)
532 :
533 : END IF
534 :
535 10180 : NULLIFY (rho_struct)
536 :
537 10180 : IF (use_virial .AND. calculate_forces) THEN
538 12 : vscale = 1.0_dp
539 : !Note: ADMMS and ADMMP stress tensor only for closed-shell calculations
540 12 : IF (admm_env%do_admms) vscale = admm_env%gsi(1)**(2.0_dp/3.0_dp)
541 12 : IF (admm_env%do_admmp) vscale = admm_env%gsi(1)**2
542 156 : virial%pv_exc = virial%pv_exc - vscale*virial%pv_xc
543 156 : virial%pv_virial = virial%pv_virial - vscale*virial%pv_xc
544 : ! virial%pv_xc will be zeroed in the xc routines
545 : END IF
546 10180 : xc_section => admm_env%xc_section_primary
547 : ELSE
548 86271 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
549 : END IF
550 :
551 96451 : IF (gapw_xc) THEN
552 2534 : CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
553 : ELSE
554 93917 : CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
555 : END IF
556 :
557 : ! zmp
558 96451 : IF (dft_control%apply_external_density .OR. dft_control%apply_external_vxc) THEN
559 0 : energy%exc = 0.0_dp
560 0 : CALL calculate_zmp_potential(qs_env, v_rspace_new, rho, exc=energy%exc)
561 : ELSE
562 : ! Embedding potential
563 96451 : IF (dft_control%apply_embed_pot) THEN
564 868 : NULLIFY (v_rspace_embed)
565 868 : energy%embed_corr = 0.0_dp
566 : CALL get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, &
567 868 : energy%embed_corr, just_energy)
568 : END IF
569 : ! Everything else
570 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
571 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=energy%exc, &
572 : edisp=edisp, dispersion_env=qs_env%dispersion_env, &
573 96451 : just_energy=just_energy_xc)
574 96451 : IF (edisp /= 0.0_dp) energy%dispersion = edisp
575 96451 : IF (qs_env%requires_matrix_vxc .AND. ASSOCIATED(v_rspace_new)) THEN
576 0 : CALL compute_matrix_vxc(qs_env=qs_env, v_rspace=v_rspace_new, matrix_vxc=matrix_vxc)
577 0 : CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
578 : END IF
579 :
580 96451 : IF (gapw .OR. gapw_xc) THEN
581 15506 : CALL calculate_vxc_atom(qs_env, just_energy_xc, energy%exc1, xc_section_external=xc_section)
582 : ! test for not implemented (bug) option
583 15506 : IF (use_virial .AND. calculate_forces) THEN
584 26 : IF (ASSOCIATED(v_tau_rspace)) THEN
585 0 : CPABORT("MGGA STRESS with GAPW/GAPW_XC not implemneted")
586 : END IF
587 : END IF
588 : END IF
589 :
590 : END IF
591 :
592 96451 : NULLIFY (rho_struct)
593 96451 : IF (use_virial .AND. calculate_forces) THEN
594 4862 : virial%pv_exc = virial%pv_exc - virial%pv_xc
595 4862 : virial%pv_virial = virial%pv_virial - virial%pv_xc
596 : END IF
597 :
598 : ! *** Add Hartree-Fock contribution if required ***
599 96451 : IF (do_hfx) THEN
600 23518 : IF (dokp) THEN
601 214 : CALL hfx_ks_matrix_kp(qs_env, ks_matrix, energy, calculate_forces)
602 : ELSE
603 : CALL hfx_ks_matrix(qs_env, ks_matrix, rho, energy, calculate_forces, &
604 23304 : just_energy, v_rspace_new, v_tau_rspace)
605 : END IF
606 : !! Adiabatic rescaling only if do_hfx; right?????
607 : END IF !do_hfx
608 :
609 96451 : IF (do_ppl .AND. calculate_forces) THEN
610 12 : CPASSERT(.NOT. gapw)
611 26 : DO ispin = 1, nspins
612 26 : CALL integrate_ppl_rspace(rho_r(ispin), qs_env)
613 : END DO
614 : END IF
615 :
616 96451 : IF (ASSOCIATED(rho_nlcc) .AND. calculate_forces) THEN
617 60 : DO ispin = 1, nspins
618 30 : CALL integrate_rho_nlcc(v_rspace_new(ispin), qs_env)
619 60 : IF (dft_control%do_admm) CALL integrate_rho_nlcc(v_rspace_new_aux_fit(ispin), qs_env)
620 : END DO
621 : END IF
622 :
623 : ! calculate KG correction
624 96451 : IF (dft_control%qs_control%do_kg .AND. just_energy) THEN
625 :
626 12 : CPASSERT(.NOT. (gapw .OR. gapw_xc))
627 12 : CPASSERT(nimages == 1)
628 12 : ksmat => ks_matrix(:, 1)
629 12 : CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces, do_kernel=.FALSE.)
630 :
631 : ! subtract kg corr from the total energy
632 12 : energy%exc = energy%exc - ekin_mol
633 :
634 : END IF
635 :
636 : ! *** Single atom contributions ***
637 96451 : IF (.NOT. just_energy) THEN
638 89468 : IF (calculate_forces) THEN
639 : ! Getting nuclear force contribution from the core charge density
640 5335 : IF ((poisson_env%parameters%solver .EQ. pw_poisson_implicit) .AND. &
641 : (poisson_env%parameters%dielectric_params%dielec_core_correction)) THEN
642 30 : BLOCK
643 : TYPE(pw_r3d_rs_type) :: v_minus_veps
644 30 : CALL auxbas_pw_pool%create_pw(v_minus_veps)
645 30 : CALL pw_copy(v_hartree_rspace, v_minus_veps)
646 30 : CALL pw_axpy(poisson_env%implicit_env%v_eps, v_minus_veps, -v_hartree_rspace%pw_grid%dvol)
647 30 : CALL integrate_v_core_rspace(v_minus_veps, qs_env)
648 60 : CALL auxbas_pw_pool%give_back_pw(v_minus_veps)
649 : END BLOCK
650 : ELSE
651 5305 : CALL integrate_v_core_rspace(v_hartree_rspace, qs_env)
652 : END IF
653 : END IF
654 :
655 89468 : IF (.NOT. do_hfx) THEN
656 : ! Initialize the Kohn-Sham matrix with the core Hamiltonian matrix
657 : ! (sets ks sparsity equal to matrix_h sparsity)
658 148335 : DO ispin = 1, nspins
659 321006 : DO img = 1, nimages
660 172671 : CALL dbcsr_get_info(ks_matrix(ispin, img)%matrix, name=name) ! keep the name
661 253288 : CALL dbcsr_copy(ks_matrix(ispin, img)%matrix, matrix_h(1, img)%matrix, name=name)
662 : END DO
663 : END DO
664 : ! imaginary part if required
665 67718 : IF (qs_env%run_rtp) THEN
666 2002 : IF (dft_control%rtp_control%velocity_gauge) THEN
667 150 : CPASSERT(ASSOCIATED(matrix_h_im))
668 150 : CPASSERT(ASSOCIATED(ks_matrix_im))
669 300 : DO ispin = 1, nspins
670 450 : DO img = 1, nimages
671 150 : CALL dbcsr_get_info(ks_matrix_im(ispin, img)%matrix, name=name) ! keep the name
672 300 : CALL dbcsr_copy(ks_matrix_im(ispin, img)%matrix, matrix_h_im(1, img)%matrix, name=name)
673 : END DO
674 : END DO
675 : END IF
676 : END IF
677 : END IF
678 :
679 89468 : IF (use_virial .AND. calculate_forces) THEN
680 4862 : pv_loc = virial%pv_virial
681 : END IF
682 : ! sum up potentials and integrate
683 : ! Pointing my_rho to the density matrix rho_ao
684 89468 : my_rho => rho_ao
685 :
686 : CALL sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, &
687 : v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, &
688 : v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, &
689 89468 : cdft_control, calculate_forces)
690 :
691 89468 : IF (use_virial .AND. calculate_forces) THEN
692 4862 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
693 : END IF
694 89468 : IF (dft_control%qs_control%do_kg) THEN
695 776 : CPASSERT(.NOT. (gapw .OR. gapw_xc))
696 776 : CPASSERT(nimages == 1)
697 776 : ksmat => ks_matrix(:, 1)
698 :
699 776 : IF (use_virial .AND. calculate_forces) THEN
700 0 : pv_loc = virial%pv_virial
701 : END IF
702 :
703 776 : CALL kg_ekin_subset(qs_env, ksmat, ekin_mol, calculate_forces, do_kernel=.FALSE.)
704 : ! subtract kg corr from the total energy
705 776 : energy%exc = energy%exc - ekin_mol
706 :
707 : ! virial corrections
708 776 : IF (use_virial .AND. calculate_forces) THEN
709 :
710 : ! Integral contribution
711 0 : virial%pv_ehartree = virial%pv_ehartree + (virial%pv_virial - pv_loc)
712 :
713 : ! GGA contribution
714 0 : virial%pv_exc = virial%pv_exc + virial%pv_xc
715 0 : virial%pv_virial = virial%pv_virial + virial%pv_xc
716 0 : virial%pv_xc = 0.0_dp
717 : END IF
718 : END IF
719 :
720 : ELSE
721 : IF (do_hfx) THEN
722 : IF (.FALSE.) THEN
723 : CPWARN("KS matrix not longer correct. Check possible problems with property calculations!")
724 : END IF
725 : END IF
726 : END IF ! .NOT. just energy
727 :
728 96451 : IF (dft_control%qs_control%ddapc_explicit_potential) THEN
729 92 : CALL auxbas_pw_pool%give_back_pw(v_spin_ddapc_rest_r)
730 92 : DEALLOCATE (v_spin_ddapc_rest_r)
731 : END IF
732 :
733 96451 : IF (calculate_forces .AND. dft_control%qs_control%cdft) THEN
734 118 : IF (.NOT. cdft_control%transfer_pot) THEN
735 212 : DO iatom = 1, SIZE(cdft_control%group)
736 114 : CALL auxbas_pw_pool%give_back_pw(cdft_control%group(iatom)%weight)
737 212 : DEALLOCATE (cdft_control%group(iatom)%weight)
738 : END DO
739 98 : IF (cdft_control%atomic_charges) THEN
740 78 : DO iatom = 1, cdft_control%natoms
741 78 : CALL auxbas_pw_pool%give_back_pw(cdft_control%charge(iatom))
742 : END DO
743 26 : DEALLOCATE (cdft_control%charge)
744 : END IF
745 98 : IF (cdft_control%type == outer_scf_becke_constraint .AND. &
746 : cdft_control%becke_control%cavity_confine) THEN
747 88 : IF (.NOT. ASSOCIATED(cdft_control%becke_control%cavity_mat)) THEN
748 64 : CALL auxbas_pw_pool%give_back_pw(cdft_control%becke_control%cavity)
749 : ELSE
750 24 : DEALLOCATE (cdft_control%becke_control%cavity_mat)
751 : END IF
752 10 : ELSE IF (cdft_control%type == outer_scf_hirshfeld_constraint) THEN
753 2 : IF (ASSOCIATED(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)) THEN
754 0 : CALL auxbas_pw_pool%give_back_pw(cdft_control%hirshfeld_control%hirshfeld_env%fnorm)
755 : END IF
756 : END IF
757 98 : IF (ASSOCIATED(cdft_control%charges_fragment)) DEALLOCATE (cdft_control%charges_fragment)
758 98 : cdft_control%save_pot = .FALSE.
759 98 : cdft_control%need_pot = .TRUE.
760 98 : cdft_control%external_control = .FALSE.
761 : END IF
762 : END IF
763 :
764 96451 : IF (dft_control%do_sccs) THEN
765 132 : CALL auxbas_pw_pool%give_back_pw(v_sccs_rspace)
766 132 : DEALLOCATE (v_sccs_rspace)
767 : END IF
768 :
769 96451 : IF (gapw) THEN
770 12972 : IF (dft_control%apply_external_potential) THEN
771 : ! Integrals of the Hartree potential with g0_soft
772 : CALL qmmm_modify_hartree_pot(v_hartree=v_hartree_rspace, &
773 42 : v_qmmm=vee, scale=-1.0_dp)
774 : END IF
775 12972 : CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces)
776 : END IF
777 :
778 96451 : IF (gapw .OR. gapw_xc) THEN
779 : ! Single atom contributions in the KS matrix ***
780 15506 : CALL update_ks_atom(qs_env, ks_matrix, rho_ao, calculate_forces)
781 15506 : IF (dft_control%do_admm) THEN
782 : !Single atom contribution to the AUX matrices
783 : !Note: also update ks_aux_fit matrix in case of rtp
784 2350 : CALL admm_update_ks_atom(qs_env, calculate_forces)
785 : END IF
786 : END IF
787 :
788 : !Calculation of Mulliken restraint, if requested
789 : CALL qs_ks_mulliken_restraint(energy, dft_control, just_energy, para_env, &
790 96451 : ks_matrix, matrix_s, rho, mulliken_order_p)
791 :
792 : ! Add DFT+U contribution, if requested
793 96451 : IF (dft_control%dft_plus_u) THEN
794 1552 : CPASSERT(nimages == 1)
795 1552 : IF (just_energy) THEN
796 588 : CALL plus_u(qs_env=qs_env)
797 : ELSE
798 964 : ksmat => ks_matrix(:, 1)
799 964 : CALL plus_u(qs_env=qs_env, matrix_h=ksmat)
800 : END IF
801 : ELSE
802 94899 : energy%dft_plus_u = 0.0_dp
803 : END IF
804 :
805 : ! At this point the ks matrix should be up to date, filter it if requested
806 212874 : DO ispin = 1, nspins
807 432425 : DO img = 1, nimages
808 : CALL dbcsr_filter(ks_matrix(ispin, img)%matrix, &
809 335974 : dft_control%qs_control%eps_filter_matrix)
810 : END DO
811 : END DO
812 :
813 : !** merge the auxiliary KS matrix and the primary one
814 96451 : IF (dft_control%do_admm_mo) THEN
815 10040 : IF (qs_env%run_rtp) THEN
816 76 : CALL rtp_admm_merge_ks_matrix(qs_env)
817 : ELSE
818 9964 : CALL admm_mo_merge_ks_matrix(qs_env)
819 : END IF
820 86411 : ELSEIF (dft_control%do_admm_dm) THEN
821 140 : CALL admm_dm_merge_ks_matrix(qs_env)
822 : END IF
823 :
824 : ! External field (nonperiodic case)
825 96451 : CALL qs_efield_local_operator(qs_env, just_energy, calculate_forces)
826 :
827 : ! Right now we can compute the orbital derivative here, as it depends currently only on the available
828 : ! Kohn-Sham matrix. This might change in the future, in which case more pieces might need to be assembled
829 : ! from this routine, notice that this part of the calculation in not linear scaling
830 : ! right now this operation is only non-trivial because of occupation numbers and the restricted keyword
831 96451 : IF (qs_env%requires_mo_derivs .AND. .NOT. just_energy .AND. .NOT. qs_env%run_rtp) THEN
832 34406 : CALL get_qs_env(qs_env, mo_derivs=mo_derivs)
833 34406 : CPASSERT(nimages == 1)
834 34406 : ksmat => ks_matrix(:, 1)
835 34406 : CALL calc_mo_derivatives(qs_env, ksmat, mo_derivs)
836 : END IF
837 :
838 : ! ADMM overlap forces
839 96451 : IF (calculate_forces .AND. dft_control%do_admm) THEN
840 262 : IF (dokp) THEN
841 24 : CALL calc_admm_ovlp_forces_kp(qs_env)
842 : ELSE
843 238 : CALL calc_admm_ovlp_forces(qs_env)
844 : END IF
845 : END IF
846 :
847 : ! deal with low spin roks
848 : CALL low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
849 96451 : calculate_forces, auxbas_pw_pool)
850 :
851 : ! deal with sic on explicit orbitals
852 : CALL sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, &
853 96451 : calculate_forces, auxbas_pw_pool)
854 :
855 : ! Periodic external field
856 96451 : CALL qs_efield_berry_phase(qs_env, just_energy, calculate_forces)
857 :
858 : ! adds s2_restraint energy and orbital derivatives
859 : CALL qs_ks_s2_restraint(dft_control, qs_env, matrix_s, &
860 96451 : energy, calculate_forces, just_energy)
861 :
862 96451 : IF (do_ppl) THEN
863 : ! update core energy for grid based local pseudopotential
864 60 : ecore_ppl = 0._dp
865 126 : DO ispin = 1, nspins
866 126 : ecore_ppl = ecore_ppl + pw_integral_ab(vppl_rspace, rho_r(ispin))
867 : END DO
868 60 : energy%core = energy%core + ecore_ppl
869 : END IF
870 :
871 96451 : IF (lrigpw) THEN
872 : ! update core energy for ppl_ri method
873 422 : CALL get_qs_env(qs_env, lri_env=lri_env, lri_density=lri_density)
874 422 : IF (lri_env%ppl_ri) THEN
875 8 : ecore_ppl = 0._dp
876 16 : DO ispin = 1, nspins
877 8 : lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
878 16 : CALL v_int_ppl_energy(qs_env, lri_v_int, ecore_ppl)
879 : END DO
880 8 : energy%core = energy%core + ecore_ppl
881 : END IF
882 : END IF
883 :
884 : ! Sum all energy terms to obtain the total energy
885 : energy%total = energy%core_overlap + energy%core_self + energy%core + energy%hartree + &
886 : energy%hartree_1c + energy%exc + energy%exc1 + energy%ex + &
887 : energy%dispersion + energy%gcp + energy%qmmm_el + energy%mulliken + &
888 : SUM(energy%ddapc_restraint) + energy%s2_restraint + &
889 : energy%dft_plus_u + energy%kTS + &
890 : energy%efield + energy%efield_core + energy%ee + &
891 : energy%ee_core + energy%exc_aux_fit + energy%image_charge + &
892 192958 : energy%sccs_pol + energy%cdft + energy%exc1_aux_fit
893 :
894 96451 : IF (dft_control%apply_embed_pot) energy%total = energy%total + energy%embed_corr
895 :
896 96451 : IF (abnormal_value(energy%total)) &
897 0 : CPABORT("KS energy is an abnormal value (NaN/Inf).")
898 :
899 : ! Print detailed energy
900 96451 : IF (my_print) THEN
901 96429 : CALL print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
902 : END IF
903 :
904 96451 : CALL timestop(handle)
905 :
906 96451 : END SUBROUTINE qs_ks_build_kohn_sham_matrix
907 :
908 : ! **************************************************************************************************
909 : !> \brief ...
910 : !> \param rho_tot_gspace ...
911 : !> \param qs_env ...
912 : !> \param rho ...
913 : !> \param skip_nuclear_density ...
914 : ! **************************************************************************************************
915 99759 : SUBROUTINE calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, skip_nuclear_density)
916 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rho_tot_gspace
917 : TYPE(qs_environment_type), POINTER :: qs_env
918 : TYPE(qs_rho_type), POINTER :: rho
919 : LOGICAL, INTENT(IN), OPTIONAL :: skip_nuclear_density
920 :
921 : INTEGER :: ispin
922 : LOGICAL :: my_skip
923 : TYPE(dft_control_type), POINTER :: dft_control
924 99759 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
925 : TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
926 : TYPE(qs_charges_type), POINTER :: qs_charges
927 :
928 99759 : my_skip = .FALSE.
929 918 : IF (PRESENT(skip_nuclear_density)) my_skip = skip_nuclear_density
930 :
931 99759 : CALL qs_rho_get(rho, rho_g=rho_g)
932 99759 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
933 :
934 99759 : IF (.NOT. my_skip) THEN
935 98851 : NULLIFY (rho_core)
936 98851 : CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
937 98851 : IF (dft_control%qs_control%gapw) THEN
938 13044 : NULLIFY (rho0_s_gs)
939 13044 : CALL get_qs_env(qs_env=qs_env, rho0_s_gs=rho0_s_gs)
940 13044 : CPASSERT(ASSOCIATED(rho0_s_gs))
941 13044 : CALL pw_copy(rho0_s_gs, rho_tot_gspace)
942 13044 : IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
943 1226 : CALL pw_axpy(rho_core, rho_tot_gspace)
944 : END IF
945 : ELSE
946 85807 : CALL pw_copy(rho_core, rho_tot_gspace)
947 : END IF
948 217960 : DO ispin = 1, dft_control%nspins
949 217960 : CALL pw_axpy(rho_g(ispin), rho_tot_gspace)
950 : END DO
951 98851 : CALL get_qs_env(qs_env=qs_env, qs_charges=qs_charges)
952 98851 : qs_charges%total_rho_gspace = pw_integrate_function(rho_tot_gspace, isign=-1)
953 : ELSE
954 1820 : DO ispin = 1, dft_control%nspins
955 1820 : CALL pw_axpy(rho_g(ispin), rho_tot_gspace)
956 : END DO
957 : END IF
958 :
959 99759 : END SUBROUTINE calc_rho_tot_gspace
960 :
961 : ! **************************************************************************************************
962 : !> \brief compute MO derivatives
963 : !> \param qs_env the qs_env to update
964 : !> \param ks_matrix ...
965 : !> \param mo_derivs ...
966 : !> \par History
967 : !> 01.2014 created, transferred from qs_ks_build_kohn_sham_matrix in
968 : !> separate subroutine
969 : !> \author Dorothea Golze
970 : ! **************************************************************************************************
971 34406 : SUBROUTINE calc_mo_derivatives(qs_env, ks_matrix, mo_derivs)
972 : TYPE(qs_environment_type), POINTER :: qs_env
973 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, mo_derivs
974 :
975 : INTEGER :: ispin
976 : LOGICAL :: uniform_occupation
977 34406 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
978 : TYPE(cp_fm_type), POINTER :: mo_coeff
979 : TYPE(dbcsr_type) :: mo_derivs2_tmp1, mo_derivs2_tmp2
980 : TYPE(dbcsr_type), POINTER :: mo_coeff_b
981 : TYPE(dft_control_type), POINTER :: dft_control
982 34406 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
983 :
984 34406 : NULLIFY (dft_control, mo_array, mo_coeff, mo_coeff_b, occupation_numbers)
985 :
986 : CALL get_qs_env(qs_env, &
987 : dft_control=dft_control, &
988 34406 : mos=mo_array)
989 :
990 75355 : DO ispin = 1, SIZE(mo_derivs)
991 :
992 : CALL get_mo_set(mo_set=mo_array(ispin), mo_coeff=mo_coeff, &
993 40949 : mo_coeff_b=mo_coeff_b, occupation_numbers=occupation_numbers)
994 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(ispin)%matrix, mo_coeff_b, &
995 40949 : 0.0_dp, mo_derivs(ispin)%matrix)
996 :
997 75355 : IF (dft_control%restricted) THEN
998 : ! only the first mo_set are actual variables, but we still need both
999 568 : CPASSERT(ispin == 1)
1000 568 : CPASSERT(SIZE(mo_array) == 2)
1001 : ! use a temporary array with the same size as the first spin for the second spin
1002 :
1003 : ! uniform_occupation is needed for this case, otherwise we can no
1004 : ! reconstruct things in ot, since we irreversibly sum
1005 568 : CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
1006 568 : CPASSERT(uniform_occupation)
1007 568 : CALL get_mo_set(mo_set=mo_array(2), uniform_occupation=uniform_occupation)
1008 568 : CPASSERT(uniform_occupation)
1009 :
1010 : ! The beta-spin might have fewer orbitals than alpa-spin...
1011 : ! create tempoary matrices with beta_nmo columns
1012 568 : CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff_b)
1013 568 : CALL dbcsr_create(mo_derivs2_tmp1, template=mo_coeff_b)
1014 :
1015 : ! calculate beta derivatives
1016 568 : CALL dbcsr_multiply('n', 'n', 1.0_dp, ks_matrix(2)%matrix, mo_coeff_b, 0.0_dp, mo_derivs2_tmp1)
1017 :
1018 : ! create larger matrix with alpha_nmo columns
1019 568 : CALL dbcsr_create(mo_derivs2_tmp2, template=mo_derivs(1)%matrix)
1020 568 : CALL dbcsr_set(mo_derivs2_tmp2, 0.0_dp)
1021 :
1022 : ! copy into larger matrix, fills the first beta_nmo columns
1023 : CALL dbcsr_copy_columns_hack(mo_derivs2_tmp2, mo_derivs2_tmp1, &
1024 : mo_array(2)%nmo, 1, 1, &
1025 : para_env=mo_array(1)%mo_coeff%matrix_struct%para_env, &
1026 568 : blacs_env=mo_array(1)%mo_coeff%matrix_struct%context)
1027 :
1028 : ! add beta contribution to alpa mo_derivs
1029 568 : CALL dbcsr_add(mo_derivs(1)%matrix, mo_derivs2_tmp2, 1.0_dp, 1.0_dp)
1030 568 : CALL dbcsr_release(mo_derivs2_tmp1)
1031 568 : CALL dbcsr_release(mo_derivs2_tmp2)
1032 : END IF
1033 : END DO
1034 :
1035 34406 : IF (dft_control%do_admm_mo) THEN
1036 5006 : CALL calc_admm_mo_derivatives(qs_env, mo_derivs)
1037 : END IF
1038 :
1039 34406 : END SUBROUTINE calc_mo_derivatives
1040 :
1041 : ! **************************************************************************************************
1042 : !> \brief updates the Kohn Sham matrix of the given qs_env (facility method)
1043 : !> \param qs_env the qs_env to update
1044 : !> \param calculate_forces if true calculate the quantities needed
1045 : !> to calculate the forces. Defaults to false.
1046 : !> \param just_energy if true updates the energies but not the
1047 : !> ks matrix. Defaults to false
1048 : !> \param print_active ...
1049 : !> \par History
1050 : !> 4.2002 created [fawzi]
1051 : !> 8.2014 kpoints [JGH]
1052 : !> 10.2014 refractored [Ole Schuett]
1053 : !> \author Fawzi Mohamed
1054 : ! **************************************************************************************************
1055 182195 : SUBROUTINE qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, &
1056 : print_active)
1057 : TYPE(qs_environment_type), POINTER :: qs_env
1058 : LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces, just_energy, &
1059 : print_active
1060 :
1061 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_ks_update_qs_env'
1062 :
1063 : INTEGER :: handle, unit_nr
1064 : LOGICAL :: c_forces, do_rebuild, energy_only, &
1065 : forces_up_to_date, potential_changed, &
1066 : rho_changed, s_mstruct_changed
1067 : TYPE(qs_ks_env_type), POINTER :: ks_env
1068 :
1069 182195 : NULLIFY (ks_env)
1070 182195 : unit_nr = cp_logger_get_default_io_unit()
1071 :
1072 182195 : c_forces = .FALSE.
1073 182195 : energy_only = .FALSE.
1074 182195 : IF (PRESENT(just_energy)) energy_only = just_energy
1075 182195 : IF (PRESENT(calculate_forces)) c_forces = calculate_forces
1076 :
1077 182195 : IF (c_forces) THEN
1078 9485 : CALL timeset(routineN//'_forces', handle)
1079 : ELSE
1080 172710 : CALL timeset(routineN, handle)
1081 : END IF
1082 :
1083 182195 : CPASSERT(ASSOCIATED(qs_env))
1084 :
1085 : CALL get_qs_env(qs_env, &
1086 : ks_env=ks_env, &
1087 : rho_changed=rho_changed, &
1088 : s_mstruct_changed=s_mstruct_changed, &
1089 : potential_changed=potential_changed, &
1090 182195 : forces_up_to_date=forces_up_to_date)
1091 :
1092 182195 : do_rebuild = .FALSE.
1093 182195 : do_rebuild = do_rebuild .OR. rho_changed
1094 7756 : do_rebuild = do_rebuild .OR. s_mstruct_changed
1095 7748 : do_rebuild = do_rebuild .OR. potential_changed
1096 7748 : do_rebuild = do_rebuild .OR. (c_forces .AND. .NOT. forces_up_to_date)
1097 :
1098 : IF (do_rebuild) THEN
1099 174797 : CALL evaluate_core_matrix_traces(qs_env)
1100 :
1101 : ! the ks matrix will be rebuilt so this is fine now
1102 174797 : CALL set_ks_env(ks_env, potential_changed=.FALSE.)
1103 :
1104 : CALL rebuild_ks_matrix(qs_env, &
1105 : calculate_forces=c_forces, &
1106 : just_energy=energy_only, &
1107 174797 : print_active=print_active)
1108 :
1109 174797 : IF (.NOT. energy_only) THEN
1110 : CALL set_ks_env(ks_env, &
1111 : rho_changed=.FALSE., &
1112 : s_mstruct_changed=.FALSE., &
1113 318683 : forces_up_to_date=forces_up_to_date .OR. c_forces)
1114 : END IF
1115 : END IF
1116 :
1117 182195 : CALL timestop(handle)
1118 :
1119 182195 : END SUBROUTINE qs_ks_update_qs_env
1120 :
1121 : ! **************************************************************************************************
1122 : !> \brief Calculates the traces of the core matrices and the density matrix.
1123 : !> \param qs_env ...
1124 : !> \author Ole Schuett
1125 : ! **************************************************************************************************
1126 174797 : SUBROUTINE evaluate_core_matrix_traces(qs_env)
1127 : TYPE(qs_environment_type), POINTER :: qs_env
1128 :
1129 : CHARACTER(LEN=*), PARAMETER :: routineN = 'evaluate_core_matrix_traces'
1130 :
1131 : INTEGER :: handle
1132 : REAL(KIND=dp) :: energy_core_im
1133 174797 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrixkp_h, matrixkp_t, rho_ao_kp
1134 : TYPE(dft_control_type), POINTER :: dft_control
1135 : TYPE(qs_energy_type), POINTER :: energy
1136 : TYPE(qs_rho_type), POINTER :: rho
1137 :
1138 174797 : CALL timeset(routineN, handle)
1139 174797 : NULLIFY (energy, rho, dft_control, rho_ao_kp, matrixkp_t, matrixkp_h)
1140 :
1141 : CALL get_qs_env(qs_env, &
1142 : rho=rho, &
1143 : energy=energy, &
1144 : dft_control=dft_control, &
1145 : kinetic_kp=matrixkp_t, &
1146 174797 : matrix_h_kp=matrixkp_h)
1147 :
1148 174797 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1149 :
1150 174797 : CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy%core, dft_control%nspins)
1151 :
1152 : ! Add the imaginary part in the RTP case
1153 174797 : IF (qs_env%run_rtp) THEN
1154 3178 : IF (dft_control%rtp_control%velocity_gauge) THEN
1155 150 : CALL get_qs_env(qs_env, matrix_h_im_kp=matrixkp_h)
1156 150 : CALL qs_rho_get(rho, rho_ao_im_kp=rho_ao_kp)
1157 150 : CALL calculate_ptrace(matrixkp_h, rho_ao_kp, energy_core_im, dft_control%nspins)
1158 150 : energy%core = energy%core - energy_core_im
1159 : END IF
1160 : END IF
1161 :
1162 : ! kinetic energy
1163 174797 : IF (ASSOCIATED(matrixkp_t)) &
1164 96323 : CALL calculate_ptrace(matrixkp_t, rho_ao_kp, energy%kinetic, dft_control%nspins)
1165 :
1166 174797 : CALL timestop(handle)
1167 174797 : END SUBROUTINE evaluate_core_matrix_traces
1168 :
1169 : ! **************************************************************************************************
1170 : !> \brief Constructs a new Khon-Sham matrix
1171 : !> \param qs_env ...
1172 : !> \param calculate_forces ...
1173 : !> \param just_energy ...
1174 : !> \param print_active ...
1175 : !> \author Ole Schuett
1176 : ! **************************************************************************************************
1177 174797 : SUBROUTINE rebuild_ks_matrix(qs_env, calculate_forces, just_energy, print_active)
1178 : TYPE(qs_environment_type), POINTER :: qs_env
1179 : LOGICAL, INTENT(IN) :: calculate_forces, just_energy
1180 : LOGICAL, INTENT(IN), OPTIONAL :: print_active
1181 :
1182 : CHARACTER(LEN=*), PARAMETER :: routineN = 'rebuild_ks_matrix'
1183 :
1184 : INTEGER :: handle
1185 : TYPE(dft_control_type), POINTER :: dft_control
1186 :
1187 174797 : CALL timeset(routineN, handle)
1188 174797 : NULLIFY (dft_control)
1189 :
1190 174797 : CALL get_qs_env(qs_env, dft_control=dft_control)
1191 :
1192 174797 : IF (dft_control%qs_control%semi_empirical) THEN
1193 : CALL build_se_fock_matrix(qs_env, &
1194 : calculate_forces=calculate_forces, &
1195 41070 : just_energy=just_energy)
1196 :
1197 133727 : ELSEIF (dft_control%qs_control%dftb) THEN
1198 : CALL build_dftb_ks_matrix(qs_env, &
1199 : calculate_forces=calculate_forces, &
1200 14400 : just_energy=just_energy)
1201 :
1202 119327 : ELSEIF (dft_control%qs_control%xtb) THEN
1203 : CALL build_xtb_ks_matrix(qs_env, &
1204 : calculate_forces=calculate_forces, &
1205 23004 : just_energy=just_energy)
1206 :
1207 : ELSE
1208 : CALL qs_ks_build_kohn_sham_matrix(qs_env, &
1209 : calculate_forces=calculate_forces, &
1210 : just_energy=just_energy, &
1211 96323 : print_active=print_active)
1212 : END IF
1213 :
1214 174797 : CALL timestop(handle)
1215 :
1216 174797 : END SUBROUTINE rebuild_ks_matrix
1217 :
1218 : ! **************************************************************************************************
1219 : !> \brief Allocate ks_matrix if necessary, take current overlap matrix as template
1220 : !> \param qs_env ...
1221 : !> \param is_complex ...
1222 : !> \par History
1223 : !> refactoring 04.03.2011 [MI]
1224 : !> \author
1225 : ! **************************************************************************************************
1226 :
1227 19552 : SUBROUTINE qs_ks_allocate_basics(qs_env, is_complex)
1228 : TYPE(qs_environment_type), POINTER :: qs_env
1229 : LOGICAL, INTENT(in) :: is_complex
1230 :
1231 : CHARACTER(LEN=default_string_length) :: headline
1232 : INTEGER :: ic, ispin, nimages, nspins
1233 : LOGICAL :: do_kpoints
1234 19552 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s_kp, matrixkp_im_ks, matrixkp_ks
1235 : TYPE(dbcsr_type), POINTER :: refmatrix
1236 : TYPE(dft_control_type), POINTER :: dft_control
1237 : TYPE(kpoint_type), POINTER :: kpoints
1238 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1239 19552 : POINTER :: sab_orb
1240 : TYPE(qs_ks_env_type), POINTER :: ks_env
1241 :
1242 19552 : NULLIFY (dft_control, ks_env, matrix_s_kp, sab_orb, matrixkp_ks, refmatrix, matrixkp_im_ks, kpoints)
1243 :
1244 : CALL get_qs_env(qs_env, &
1245 : dft_control=dft_control, &
1246 : matrix_s_kp=matrix_s_kp, &
1247 : ks_env=ks_env, &
1248 : kpoints=kpoints, &
1249 : do_kpoints=do_kpoints, &
1250 : matrix_ks_kp=matrixkp_ks, &
1251 19552 : matrix_ks_im_kp=matrixkp_im_ks)
1252 :
1253 19552 : IF (do_kpoints) THEN
1254 764 : CALL get_kpoint_info(kpoints, sab_nl=sab_orb)
1255 : ELSE
1256 18788 : CALL get_qs_env(qs_env, sab_orb=sab_orb)
1257 : END IF
1258 :
1259 19552 : nspins = dft_control%nspins
1260 19552 : nimages = dft_control%nimages
1261 :
1262 19552 : IF (.NOT. ASSOCIATED(matrixkp_ks)) THEN
1263 19512 : CALL dbcsr_allocate_matrix_set(matrixkp_ks, nspins, nimages)
1264 19512 : refmatrix => matrix_s_kp(1, 1)%matrix
1265 41508 : DO ispin = 1, nspins
1266 132692 : DO ic = 1, nimages
1267 91184 : IF (nspins > 1) THEN
1268 24680 : IF (ispin == 1) THEN
1269 12340 : headline = "KOHN-SHAM MATRIX FOR ALPHA SPIN"
1270 : ELSE
1271 12340 : headline = "KOHN-SHAM MATRIX FOR BETA SPIN"
1272 : END IF
1273 : ELSE
1274 66504 : headline = "KOHN-SHAM MATRIX"
1275 : END IF
1276 91184 : ALLOCATE (matrixkp_ks(ispin, ic)%matrix)
1277 : CALL dbcsr_create(matrix=matrixkp_ks(ispin, ic)%matrix, template=refmatrix, &
1278 91184 : name=TRIM(headline), matrix_type=dbcsr_type_symmetric, nze=0)
1279 91184 : CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_ks(ispin, ic)%matrix, sab_orb)
1280 113180 : CALL dbcsr_set(matrixkp_ks(ispin, ic)%matrix, 0.0_dp)
1281 : END DO
1282 : END DO
1283 19512 : CALL set_ks_env(ks_env, matrix_ks_kp=matrixkp_ks)
1284 : END IF
1285 :
1286 19552 : IF (is_complex) THEN
1287 138 : IF (.NOT. ASSOCIATED(matrixkp_im_ks)) THEN
1288 138 : CPASSERT(nspins .EQ. SIZE(matrixkp_ks, 1))
1289 138 : CPASSERT(nimages .EQ. SIZE(matrixkp_ks, 2))
1290 138 : CALL dbcsr_allocate_matrix_set(matrixkp_im_ks, nspins, nimages)
1291 288 : DO ispin = 1, nspins
1292 438 : DO ic = 1, nimages
1293 150 : IF (nspins > 1) THEN
1294 24 : IF (ispin == 1) THEN
1295 12 : headline = "IMAGINARY KOHN-SHAM MATRIX FOR ALPHA SPIN"
1296 : ELSE
1297 12 : headline = "IMAGINARY KOHN-SHAM MATRIX FOR BETA SPIN"
1298 : END IF
1299 : ELSE
1300 126 : headline = "IMAGINARY KOHN-SHAM MATRIX"
1301 : END IF
1302 150 : ALLOCATE (matrixkp_im_ks(ispin, ic)%matrix)
1303 150 : refmatrix => matrixkp_ks(ispin, ic)%matrix ! base on real part, but anti-symmetric
1304 : CALL dbcsr_create(matrix=matrixkp_im_ks(ispin, ic)%matrix, template=refmatrix, &
1305 150 : name=TRIM(headline), matrix_type=dbcsr_type_antisymmetric, nze=0)
1306 150 : CALL cp_dbcsr_alloc_block_from_nbl(matrixkp_im_ks(ispin, ic)%matrix, sab_orb)
1307 300 : CALL dbcsr_set(matrixkp_im_ks(ispin, ic)%matrix, 0.0_dp)
1308 : END DO
1309 : END DO
1310 138 : CALL set_ks_env(ks_env, matrix_ks_im_kp=matrixkp_im_ks)
1311 : END IF
1312 : END IF
1313 :
1314 19552 : END SUBROUTINE qs_ks_allocate_basics
1315 :
1316 : END MODULE qs_ks_methods
|
