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