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