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