(git:15c1bfc)
Loading...
Searching...
No Matches
qs_ks_utils.F
Go to the documentation of this file.
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!> \par History
12!> 05.2002 moved from qs_scf (see there the history) [fawzi]
13!> JGH [30.08.02] multi-grid arrays independent from density and potential
14!> 10.2002 introduced pools, uses updated rho as input,
15!> removed most temporary variables, renamed may vars,
16!> began conversion to LSD [fawzi]
17!> 10.2004 moved calculate_w_matrix here [Joost VandeVondele]
18!> introduced energy derivative wrt MOs [Joost VandeVondele]
19!> \author Fawzi Mohamed
20! **************************************************************************************************
21
23 USE admm_types, ONLY: admm_type,&
26 USE cell_types, ONLY: cell_type
28 USE cp_dbcsr_api, ONLY: &
31 USE cp_dbcsr_contrib, ONLY: dbcsr_dot,&
43 USE cp_fm_types, ONLY: cp_fm_create,&
52 USE cp_output_handling, ONLY: cp_p_file,&
58 USE hfx_types, ONLY: hfx_type
59 USE input_constants, ONLY: &
68 USE kinds, ONLY: default_string_length,&
69 dp
70 USE kpoint_types, ONLY: get_kpoint_info,&
81 USE ps_implicit_types, ONLY: mixed_bc,&
85 USE pw_env_types, ONLY: pw_env_get,&
87 USE pw_methods, ONLY: pw_axpy,&
88 pw_copy,&
91 pw_scale,&
98 USE pw_types, ONLY: pw_c1d_gs_type,&
107 USE qs_integrate_potential, ONLY: integrate_v_rspace,&
108 integrate_v_rspace_diagonal,&
109 integrate_v_rspace_one_center
110 USE qs_kind_types, ONLY: get_qs_kind_set,&
113 USE qs_ks_types, ONLY: get_ks_env,&
115 USE qs_mo_types, ONLY: get_mo_set,&
117 USE qs_rho_types, ONLY: qs_rho_get,&
120 USE virial_types, ONLY: virial_type
121 USE xc, ONLY: xc_exc_calc,&
123#include "./base/base_uses.f90"
124
125 IMPLICIT NONE
126
127 PRIVATE
128
129 LOGICAL, PARAMETER :: debug_this_module = .true.
130 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_utils'
131
135
136CONTAINS
137
138! **************************************************************************************************
139!> \brief do ROKS calculations yielding low spin states
140!> \param energy ...
141!> \param qs_env ...
142!> \param dft_control ...
143!> \param do_hfx ...
144!> \param just_energy ...
145!> \param calculate_forces ...
146!> \param auxbas_pw_pool ...
147! **************************************************************************************************
148 SUBROUTINE low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
149 calculate_forces, auxbas_pw_pool)
150
151 TYPE(qs_energy_type), POINTER :: energy
152 TYPE(qs_environment_type), POINTER :: qs_env
153 TYPE(dft_control_type), POINTER :: dft_control
154 LOGICAL, INTENT(IN) :: do_hfx, just_energy, calculate_forces
155 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
156
157 CHARACTER(*), PARAMETER :: routinen = 'low_spin_roks'
158
159 INTEGER :: handle, irep, ispin, iterm, k, k_alpha, &
160 k_beta, n_rep, nelectron, nspin, nterms
161 INTEGER, DIMENSION(:), POINTER :: ivec
162 INTEGER, DIMENSION(:, :, :), POINTER :: occupations
163 LOGICAL :: compute_virial, in_range, &
164 uniform_occupation
165 REAL(kind=dp) :: ehfx, exc
166 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
167 REAL(kind=dp), DIMENSION(:), POINTER :: energy_scaling, rvec, scaling
168 TYPE(cell_type), POINTER :: cell
169 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_hfx, matrix_p, mdummy, &
170 mo_derivs, rho_ao
171 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p2
172 TYPE(dbcsr_type), POINTER :: dbcsr_deriv, fm_deriv, fm_scaled, &
173 mo_coeff
174 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
175 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
176 TYPE(mp_para_env_type), POINTER :: para_env
177 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
178 TYPE(pw_env_type), POINTER :: pw_env
179 TYPE(pw_pool_type), POINTER :: xc_pw_pool
180 TYPE(pw_r3d_rs_type) :: work_v_rspace
181 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc, vxc_tau
182 TYPE(qs_ks_env_type), POINTER :: ks_env
183 TYPE(qs_rho_type), POINTER :: rho
184 TYPE(section_vals_type), POINTER :: hfx_section, input, &
185 low_spin_roks_section, xc_section
186 TYPE(virial_type), POINTER :: virial
187
188 IF (.NOT. dft_control%low_spin_roks) RETURN
189
190 CALL timeset(routinen, handle)
191
192 NULLIFY (ks_env, rho_ao)
193
194 ! Test for not compatible options
195 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
196 CALL cp_abort(__location__, "GAPW/GAPW_XC are not compatible with low spin ROKS method.")
197 END IF
198 IF (dft_control%do_admm) THEN
199 CALL cp_abort(__location__, "ADMM not compatible with low spin ROKS method.")
200 END IF
201 IF (dft_control%do_admm) THEN
202 IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
203 CALL cp_abort(__location__, "ADMM with XC correction functional "// &
204 "not compatible with low spin ROKS method.")
205 END IF
206 END IF
207 IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
208 dft_control%qs_control%xtb) THEN
209 CALL cp_abort(__location__, "SE/xTB/DFTB are not compatible with low spin ROKS method.")
210 END IF
211
212 CALL get_qs_env(qs_env, &
213 ks_env=ks_env, &
214 mo_derivs=mo_derivs, &
215 mos=mo_array, &
216 rho=rho, &
217 pw_env=pw_env, &
218 input=input, &
219 cell=cell, &
220 virial=virial)
221
222 CALL qs_rho_get(rho, rho_ao=rho_ao)
223
224 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
225 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
226 hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
227
228 ! some assumptions need to be checked
229 ! we have two spins
230 cpassert(SIZE(mo_array, 1) == 2)
231 nspin = 2
232 ! we want uniform occupations
233 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
234 cpassert(uniform_occupation)
235 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation)
236 cpassert(uniform_occupation)
237 IF (do_hfx .AND. calculate_forces .AND. compute_virial) THEN
238 CALL cp_abort(__location__, "ROKS virial with HFX not available.")
239 END IF
240
241 NULLIFY (dbcsr_deriv)
242 CALL dbcsr_init_p(dbcsr_deriv)
243 CALL dbcsr_copy(dbcsr_deriv, mo_derivs(1)%matrix)
244 CALL dbcsr_set(dbcsr_deriv, 0.0_dp)
245
246 ! basic info
247 CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
248 CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_alpha)
249 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff)
250 CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_beta)
251
252 ! read the input
253 low_spin_roks_section => section_vals_get_subs_vals(input, "DFT%LOW_SPIN_ROKS")
254
255 CALL section_vals_val_get(low_spin_roks_section, "ENERGY_SCALING", r_vals=rvec)
256 nterms = SIZE(rvec)
257 ALLOCATE (energy_scaling(nterms))
258 energy_scaling = rvec !? just wondering, should this add up to 1, in which case we should cpp?
259
260 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", n_rep_val=n_rep)
261 cpassert(n_rep == nterms)
262 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=1, i_vals=ivec)
263 nelectron = SIZE(ivec)
264 cpassert(nelectron == k_alpha - k_beta)
265 ALLOCATE (occupations(2, nelectron, nterms))
266 occupations = 0
267 DO iterm = 1, nterms
268 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=iterm, i_vals=ivec)
269 cpassert(nelectron == SIZE(ivec))
270 in_range = all(ivec >= 1) .AND. all(ivec <= 2)
271 cpassert(in_range)
272 DO k = 1, nelectron
273 occupations(ivec(k), k, iterm) = 1
274 END DO
275 END DO
276
277 ! set up general data structures
278 ! density matrices, kohn-sham matrices
279
280 NULLIFY (matrix_p)
281 CALL dbcsr_allocate_matrix_set(matrix_p, nspin)
282 DO ispin = 1, nspin
283 ALLOCATE (matrix_p(ispin)%matrix)
284 CALL dbcsr_copy(matrix_p(ispin)%matrix, rho_ao(1)%matrix, &
285 name="density matrix low spin roks")
286 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
287 END DO
288
289 NULLIFY (matrix_h)
290 CALL dbcsr_allocate_matrix_set(matrix_h, nspin)
291 DO ispin = 1, nspin
292 ALLOCATE (matrix_h(ispin)%matrix)
293 CALL dbcsr_copy(matrix_h(ispin)%matrix, rho_ao(1)%matrix, &
294 name="KS matrix low spin roks")
295 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
296 END DO
297
298 IF (do_hfx) THEN
299 NULLIFY (matrix_hfx)
300 CALL dbcsr_allocate_matrix_set(matrix_hfx, nspin)
301 DO ispin = 1, nspin
302 ALLOCATE (matrix_hfx(ispin)%matrix)
303 CALL dbcsr_copy(matrix_hfx(ispin)%matrix, rho_ao(1)%matrix, &
304 name="HFX matrix low spin roks")
305 END DO
306 END IF
307
308 ! grids in real and g space for rho and vxc
309 ! tau functionals are not supported
310 NULLIFY (tau, vxc_tau, vxc)
311 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
312
313 ALLOCATE (rho_r(nspin))
314 ALLOCATE (rho_g(nspin))
315 DO ispin = 1, nspin
316 CALL auxbas_pw_pool%create_pw(rho_r(ispin))
317 CALL auxbas_pw_pool%create_pw(rho_g(ispin))
318 END DO
319 CALL auxbas_pw_pool%create_pw(work_v_rspace)
320
321 ! get mo matrices needed to construct the density matrices
322 ! we will base all on the alpha spin matrix, obviously possible in ROKS
323 CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
324 NULLIFY (fm_scaled, fm_deriv)
325 CALL dbcsr_init_p(fm_scaled)
326 CALL dbcsr_init_p(fm_deriv)
327 CALL dbcsr_copy(fm_scaled, mo_coeff)
328 CALL dbcsr_copy(fm_deriv, mo_coeff)
329
330 ALLOCATE (scaling(k_alpha))
331
332 ! for each term, add it with the given scaling factor to the energy, and compute the required derivatives
333 DO iterm = 1, nterms
334
335 DO ispin = 1, nspin
336 ! compute the proper density matrices with the required occupations
337 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
338 scaling = 1.0_dp
339 scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
340 CALL dbcsr_copy(fm_scaled, mo_coeff)
341 CALL dbcsr_scale_by_vector(fm_scaled, scaling, side='right')
342 CALL dbcsr_multiply('n', 't', 1.0_dp, mo_coeff, fm_scaled, &
343 0.0_dp, matrix_p(ispin)%matrix, retain_sparsity=.true.)
344 ! compute the densities on the grid
345 CALL calculate_rho_elec(matrix_p=matrix_p(ispin)%matrix, &
346 rho=rho_r(ispin), rho_gspace=rho_g(ispin), &
347 ks_env=ks_env)
348 END DO
349
350 ! compute the exchange energies / potential if needed
351 IF (just_energy) THEN
352 exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
353 pw_pool=xc_pw_pool)
354 ELSE
355 cpassert(.NOT. compute_virial)
356 CALL xc_vxc_pw_create1(vxc_rho=vxc, rho_r=rho_r, &
357 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
358 pw_pool=xc_pw_pool, compute_virial=.false., virial_xc=virial_xc_tmp)
359 END IF
360
361 energy%exc = energy%exc + energy_scaling(iterm)*exc
362
363 IF (do_hfx) THEN
364 ! Add Hartree-Fock contribution
365 DO ispin = 1, nspin
366 CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
367 END DO
368 ehfx = energy%ex
369 CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, &
370 recalc_integrals=.false., update_energy=.true.)
371 energy%ex = ehfx + energy_scaling(iterm)*energy%ex
372 END IF
373
374 ! add the corresponding derivatives to the MO derivatives
375 IF (.NOT. just_energy) THEN
376 ! get the potential in matrix form
377 DO ispin = 1, nspin
378 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
379 ! use a work_v_rspace
380 CALL pw_axpy(vxc(ispin), work_v_rspace, energy_scaling(iterm)*vxc(ispin)%pw_grid%dvol, 0.0_dp)
381 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), &
382 qs_env=qs_env, calculate_forces=calculate_forces)
383 CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
384 END DO
385 DEALLOCATE (vxc)
386
387 IF (do_hfx) THEN
388 ! add HFX contribution
389 DO ispin = 1, nspin
390 CALL dbcsr_add(matrix_h(ispin)%matrix, matrix_hfx(ispin)%matrix, &
391 1.0_dp, energy_scaling(iterm))
392 END DO
393 IF (calculate_forces) THEN
394 CALL get_qs_env(qs_env, x_data=x_data, para_env=para_env)
395 IF (x_data(1, 1)%n_rep_hf /= 1) THEN
396 CALL cp_abort(__location__, "Multiple HFX section forces not compatible "// &
397 "with low spin ROKS method.")
398 END IF
399 IF (x_data(1, 1)%do_hfx_ri) THEN
400 CALL cp_abort(__location__, "HFX_RI forces not compatible with low spin ROKS method.")
401 ELSE
402 irep = 1
403 NULLIFY (mdummy)
404 matrix_p2(1:nspin, 1:1) => matrix_p(1:nspin)
405 CALL derivatives_four_center(qs_env, matrix_p2, mdummy, hfx_section, para_env, &
406 irep, compute_virial, &
407 adiabatic_rescale_factor=energy_scaling(iterm))
408 END IF
409 END IF
410
411 END IF
412
413 ! add this to the mo_derivs, again based on the alpha mo_coeff
414 DO ispin = 1, nspin
415 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, &
416 0.0_dp, dbcsr_deriv, last_column=k_alpha)
417
418 scaling = 1.0_dp
419 scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
420 CALL dbcsr_scale_by_vector(dbcsr_deriv, scaling, side='right')
421 CALL dbcsr_add(mo_derivs(1)%matrix, dbcsr_deriv, 1.0_dp, 1.0_dp)
422 END DO
423
424 END IF
425
426 END DO
427
428 ! release allocated memory
429 DO ispin = 1, nspin
430 CALL auxbas_pw_pool%give_back_pw(rho_r(ispin))
431 CALL auxbas_pw_pool%give_back_pw(rho_g(ispin))
432 END DO
433 DEALLOCATE (rho_r, rho_g)
434 CALL dbcsr_deallocate_matrix_set(matrix_p)
435 CALL dbcsr_deallocate_matrix_set(matrix_h)
436 IF (do_hfx) THEN
437 CALL dbcsr_deallocate_matrix_set(matrix_hfx)
438 END IF
439
440 CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
441
442 CALL dbcsr_release_p(fm_deriv)
443 CALL dbcsr_release_p(fm_scaled)
444
445 DEALLOCATE (occupations)
446 DEALLOCATE (energy_scaling)
447 DEALLOCATE (scaling)
448
449 CALL dbcsr_release_p(dbcsr_deriv)
450
451 CALL timestop(handle)
452
453 END SUBROUTINE low_spin_roks
454! **************************************************************************************************
455!> \brief do sic calculations on explicit orbitals
456!> \param energy ...
457!> \param qs_env ...
458!> \param dft_control ...
459!> \param poisson_env ...
460!> \param just_energy ...
461!> \param calculate_forces ...
462!> \param auxbas_pw_pool ...
463! **************************************************************************************************
464 SUBROUTINE sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, &
465 calculate_forces, auxbas_pw_pool)
466
467 TYPE(qs_energy_type), POINTER :: energy
468 TYPE(qs_environment_type), POINTER :: qs_env
469 TYPE(dft_control_type), POINTER :: dft_control
470 TYPE(pw_poisson_type), POINTER :: poisson_env
471 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
472 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
473
474 CHARACTER(*), PARAMETER :: routinen = 'sic_explicit_orbitals'
475
476 INTEGER :: handle, i, iorb, k_alpha, k_beta, norb
477 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: sic_orbital_list
478 LOGICAL :: compute_virial, uniform_occupation
479 REAL(kind=dp) :: ener, exc
480 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
481 TYPE(cell_type), POINTER :: cell
482 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
483 TYPE(cp_fm_type) :: matrix_hv, matrix_v
484 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_derivs_local
485 TYPE(cp_fm_type), POINTER :: mo_coeff
486 TYPE(dbcsr_p_type) :: orb_density_matrix_p, orb_h_p
487 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, rho_ao, tmp_dbcsr
488 TYPE(dbcsr_type), POINTER :: orb_density_matrix, orb_h
489 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
490 TYPE(pw_c1d_gs_type) :: work_v_gspace
491 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
492 TYPE(pw_c1d_gs_type), TARGET :: orb_rho_g, tmp_g
493 TYPE(pw_env_type), POINTER :: pw_env
494 TYPE(pw_pool_type), POINTER :: xc_pw_pool
495 TYPE(pw_r3d_rs_type) :: work_v_rspace
496 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc, vxc_tau
497 TYPE(pw_r3d_rs_type), TARGET :: orb_rho_r, tmp_r
498 TYPE(qs_ks_env_type), POINTER :: ks_env
499 TYPE(qs_rho_type), POINTER :: rho
500 TYPE(section_vals_type), POINTER :: input, xc_section
501 TYPE(virial_type), POINTER :: virial
502
503 IF (dft_control%sic_method_id .NE. sic_eo) RETURN
504
505 CALL timeset(routinen, handle)
506
507 NULLIFY (tau, vxc_tau, mo_derivs, ks_env, rho_ao)
508
509 ! generate the lists of orbitals that need sic treatment
510 CALL get_qs_env(qs_env, &
511 ks_env=ks_env, &
512 mo_derivs=mo_derivs, &
513 mos=mo_array, &
514 rho=rho, &
515 pw_env=pw_env, &
516 input=input, &
517 cell=cell, &
518 virial=virial)
519
520 CALL qs_rho_get(rho, rho_ao=rho_ao)
521
522 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
523 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
524
525 DO i = 1, SIZE(mo_array) !fm->dbcsr
526 IF (mo_array(i)%use_mo_coeff_b) THEN !fm->dbcsr
527 CALL copy_dbcsr_to_fm(mo_array(i)%mo_coeff_b, &
528 mo_array(i)%mo_coeff) !fm->dbcsr
529 END IF !fm->dbcsr
530 END DO !fm->dbcsr
531
532 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
533
534 ! we have two spins
535 cpassert(SIZE(mo_array, 1) == 2)
536 ! we want uniform occupations
537 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
538 cpassert(uniform_occupation)
539 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff, uniform_occupation=uniform_occupation)
540 cpassert(uniform_occupation)
541
542 NULLIFY (tmp_dbcsr)
543 CALL dbcsr_allocate_matrix_set(tmp_dbcsr, SIZE(mo_derivs, 1))
544 DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
545 !
546 NULLIFY (tmp_dbcsr(i)%matrix)
547 CALL dbcsr_init_p(tmp_dbcsr(i)%matrix)
548 CALL dbcsr_copy(tmp_dbcsr(i)%matrix, mo_derivs(i)%matrix)
549 CALL dbcsr_set(tmp_dbcsr(i)%matrix, 0.0_dp)
550 END DO !fm->dbcsr
551
552 k_alpha = 0; k_beta = 0
553 SELECT CASE (dft_control%sic_list_id)
554 CASE (sic_list_all)
555
556 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
557 CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
558
559 IF (SIZE(mo_array, 1) > 1) THEN
560 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
561 CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
562 END IF
563
564 norb = k_alpha + k_beta
565 ALLOCATE (sic_orbital_list(3, norb))
566
567 iorb = 0
568 DO i = 1, k_alpha
569 iorb = iorb + 1
570 sic_orbital_list(1, iorb) = 1
571 sic_orbital_list(2, iorb) = i
572 sic_orbital_list(3, iorb) = 1
573 END DO
574 DO i = 1, k_beta
575 iorb = iorb + 1
576 sic_orbital_list(1, iorb) = 2
577 sic_orbital_list(2, iorb) = i
578 IF (SIZE(mo_derivs, 1) == 1) THEN
579 sic_orbital_list(3, iorb) = 1
580 ELSE
581 sic_orbital_list(3, iorb) = 2
582 END IF
583 END DO
584
585 CASE (sic_list_unpaired)
586 ! we have two spins
587 cpassert(SIZE(mo_array, 1) == 2)
588 ! we have them restricted
589 cpassert(SIZE(mo_derivs, 1) == 1)
590 cpassert(dft_control%restricted)
591
592 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
593 CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
594
595 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
596 CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
597
598 norb = k_alpha - k_beta
599 ALLOCATE (sic_orbital_list(3, norb))
600
601 iorb = 0
602 DO i = k_beta + 1, k_alpha
603 iorb = iorb + 1
604 sic_orbital_list(1, iorb) = 1
605 sic_orbital_list(2, iorb) = i
606 ! we are guaranteed to be restricted
607 sic_orbital_list(3, iorb) = 1
608 END DO
609
610 CASE DEFAULT
611 cpabort("")
612 END SELECT
613
614 ! data needed for each of the orbs
615 CALL auxbas_pw_pool%create_pw(orb_rho_r)
616 CALL auxbas_pw_pool%create_pw(tmp_r)
617 CALL auxbas_pw_pool%create_pw(orb_rho_g)
618 CALL auxbas_pw_pool%create_pw(tmp_g)
619 CALL auxbas_pw_pool%create_pw(work_v_gspace)
620 CALL auxbas_pw_pool%create_pw(work_v_rspace)
621
622 ALLOCATE (orb_density_matrix)
623 CALL dbcsr_copy(orb_density_matrix, rho_ao(1)%matrix, &
624 name="orb_density_matrix")
625 CALL dbcsr_set(orb_density_matrix, 0.0_dp)
626 orb_density_matrix_p%matrix => orb_density_matrix
627
628 ALLOCATE (orb_h)
629 CALL dbcsr_copy(orb_h, rho_ao(1)%matrix, &
630 name="orb_density_matrix")
631 CALL dbcsr_set(orb_h, 0.0_dp)
632 orb_h_p%matrix => orb_h
633
634 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
635
636 CALL cp_fm_struct_create(fm_struct_tmp, ncol_global=1, &
637 template_fmstruct=mo_coeff%matrix_struct)
638 CALL cp_fm_create(matrix_v, fm_struct_tmp, name="matrix_v")
639 CALL cp_fm_create(matrix_hv, fm_struct_tmp, name="matrix_hv")
640 CALL cp_fm_struct_release(fm_struct_tmp)
641
642 ALLOCATE (mo_derivs_local(SIZE(mo_array, 1)))
643 DO i = 1, SIZE(mo_array, 1)
644 CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff)
645 CALL cp_fm_create(mo_derivs_local(i), mo_coeff%matrix_struct)
646 END DO
647
648 ALLOCATE (rho_r(2))
649 rho_r(1) = orb_rho_r
650 rho_r(2) = tmp_r
651 CALL pw_zero(tmp_r)
652
653 ALLOCATE (rho_g(2))
654 rho_g(1) = orb_rho_g
655 rho_g(2) = tmp_g
656 CALL pw_zero(tmp_g)
657
658 NULLIFY (vxc)
659 ! now apply to SIC correction to each selected orbital
660 DO iorb = 1, norb
661 ! extract the proper orbital from the mo_coeff
662 CALL get_mo_set(mo_set=mo_array(sic_orbital_list(1, iorb)), mo_coeff=mo_coeff)
663 CALL cp_fm_to_fm(mo_coeff, matrix_v, 1, sic_orbital_list(2, iorb), 1)
664
665 ! construct the density matrix and the corresponding density
666 CALL dbcsr_set(orb_density_matrix, 0.0_dp)
667 CALL cp_dbcsr_plus_fm_fm_t(orb_density_matrix, matrix_v=matrix_v, ncol=1, &
668 alpha=1.0_dp)
669
670 CALL calculate_rho_elec(matrix_p=orb_density_matrix, &
671 rho=orb_rho_r, rho_gspace=orb_rho_g, &
672 ks_env=ks_env)
673
674 ! compute the energy functional for this orbital and its derivative
675
676 CALL pw_poisson_solve(poisson_env, orb_rho_g, ener, work_v_gspace)
677 ! no PBC correction is done here, see "calc_v_sic_rspace" for SIC methods
678 ! with PBC aware corrections
679 energy%hartree = energy%hartree - dft_control%sic_scaling_a*ener
680 IF (.NOT. just_energy) THEN
681 CALL pw_transfer(work_v_gspace, work_v_rspace)
682 CALL pw_scale(work_v_rspace, -dft_control%sic_scaling_a*work_v_rspace%pw_grid%dvol)
683 CALL dbcsr_set(orb_h, 0.0_dp)
684 END IF
685
686 IF (just_energy) THEN
687 exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
688 pw_pool=xc_pw_pool)
689 ELSE
690 cpassert(.NOT. compute_virial)
691 CALL xc_vxc_pw_create1(vxc_rho=vxc, rho_r=rho_r, &
692 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
693 pw_pool=xc_pw_pool, compute_virial=compute_virial, virial_xc=virial_xc_tmp)
694 ! add to the existing work_v_rspace
695 CALL pw_axpy(vxc(1), work_v_rspace, -dft_control%sic_scaling_b*vxc(1)%pw_grid%dvol)
696 END IF
697 energy%exc = energy%exc - dft_control%sic_scaling_b*exc
698
699 IF (.NOT. just_energy) THEN
700 ! note, orb_h (which is being pointed to with orb_h_p) is zeroed above
701 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=orb_density_matrix_p, hmat=orb_h_p, &
702 qs_env=qs_env, calculate_forces=calculate_forces)
703
704 ! add this to the mo_derivs
705 CALL cp_dbcsr_sm_fm_multiply(orb_h, matrix_v, matrix_hv, 1)
706 ! silly trick, copy to an array of the right size and add to mo_derivs
707 CALL cp_fm_set_all(mo_derivs_local(sic_orbital_list(3, iorb)), 0.0_dp)
708 CALL cp_fm_to_fm(matrix_hv, mo_derivs_local(sic_orbital_list(3, iorb)), 1, 1, sic_orbital_list(2, iorb))
709 CALL copy_fm_to_dbcsr(mo_derivs_local(sic_orbital_list(3, iorb)), &
710 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix)
711 CALL dbcsr_add(mo_derivs(sic_orbital_list(3, iorb))%matrix, &
712 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix, 1.0_dp, 1.0_dp)
713 !
714 ! need to deallocate vxc
715 CALL xc_pw_pool%give_back_pw(vxc(1))
716 CALL xc_pw_pool%give_back_pw(vxc(2))
717 DEALLOCATE (vxc)
718
719 END IF
720
721 END DO
722
723 CALL auxbas_pw_pool%give_back_pw(orb_rho_r)
724 CALL auxbas_pw_pool%give_back_pw(tmp_r)
725 CALL auxbas_pw_pool%give_back_pw(orb_rho_g)
726 CALL auxbas_pw_pool%give_back_pw(tmp_g)
727 CALL auxbas_pw_pool%give_back_pw(work_v_gspace)
728 CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
729
730 CALL dbcsr_deallocate_matrix(orb_density_matrix)
731 CALL dbcsr_deallocate_matrix(orb_h)
732 CALL cp_fm_release(matrix_v)
733 CALL cp_fm_release(matrix_hv)
734 CALL cp_fm_release(mo_derivs_local)
735 DEALLOCATE (rho_r)
736 DEALLOCATE (rho_g)
737
738 CALL dbcsr_deallocate_matrix_set(tmp_dbcsr) !fm->dbcsr
739
740 CALL timestop(handle)
741
742 END SUBROUTINE sic_explicit_orbitals
743
744! **************************************************************************************************
745!> \brief do sic calculations on the spin density
746!> \param v_sic_rspace ...
747!> \param energy ...
748!> \param qs_env ...
749!> \param dft_control ...
750!> \param rho ...
751!> \param poisson_env ...
752!> \param just_energy ...
753!> \param calculate_forces ...
754!> \param auxbas_pw_pool ...
755! **************************************************************************************************
756 SUBROUTINE calc_v_sic_rspace(v_sic_rspace, energy, &
757 qs_env, dft_control, rho, poisson_env, just_energy, &
758 calculate_forces, auxbas_pw_pool)
759
760 TYPE(pw_r3d_rs_type), POINTER :: v_sic_rspace
761 TYPE(qs_energy_type), POINTER :: energy
762 TYPE(qs_environment_type), POINTER :: qs_env
763 TYPE(dft_control_type), POINTER :: dft_control
764 TYPE(qs_rho_type), POINTER :: rho
765 TYPE(pw_poisson_type), POINTER :: poisson_env
766 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
767 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
768
769 INTEGER :: i, nelec, nelec_a, nelec_b, nforce
770 REAL(kind=dp) :: ener, full_scaling, scaling
771 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: store_forces
772 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
773 TYPE(pw_c1d_gs_type) :: work_rho, work_v
774 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
775 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
776
777 NULLIFY (mo_array, rho_g)
778
779 IF (dft_control%sic_method_id == sic_none) RETURN
780 IF (dft_control%sic_method_id == sic_eo) RETURN
781
782 IF (dft_control%qs_control%gapw) &
783 cpabort("sic and GAPW not yet compatible")
784
785 ! OK, right now we like two spins to do sic, could be relaxed for AD
786 cpassert(dft_control%nspins == 2)
787
788 CALL auxbas_pw_pool%create_pw(work_rho)
789 CALL auxbas_pw_pool%create_pw(work_v)
790
791 CALL qs_rho_get(rho, rho_g=rho_g)
792
793 ! Hartree sic corrections
794 SELECT CASE (dft_control%sic_method_id)
796 CALL pw_copy(rho_g(1), work_rho)
797 CALL pw_axpy(rho_g(2), work_rho, alpha=-1._dp)
798 CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
799 CASE (sic_ad)
800 ! find out how many elecs we have
801 CALL get_qs_env(qs_env, mos=mo_array)
802 CALL get_mo_set(mo_set=mo_array(1), nelectron=nelec_a)
803 CALL get_mo_set(mo_set=mo_array(2), nelectron=nelec_b)
804 nelec = nelec_a + nelec_b
805 CALL pw_copy(rho_g(1), work_rho)
806 CALL pw_axpy(rho_g(2), work_rho)
807 scaling = 1.0_dp/real(nelec, kind=dp)
808 CALL pw_scale(work_rho, scaling)
809 CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
810 CASE DEFAULT
811 cpabort("Unknown sic method id")
812 END SELECT
813
814 ! Correct for DDAP charges (if any)
815 ! storing whatever force might be there from previous decoupling
816 IF (calculate_forces) THEN
817 CALL get_qs_env(qs_env=qs_env, force=force)
818 nforce = 0
819 DO i = 1, SIZE(force)
820 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
821 END DO
822 ALLOCATE (store_forces(3, nforce))
823 nforce = 0
824 DO i = 1, SIZE(force)
825 store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2)) = force(i)%ch_pulay(:, :)
826 force(i)%ch_pulay(:, :) = 0.0_dp
827 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
828 END DO
829 END IF
830
831 CALL cp_ddapc_apply_cd(qs_env, &
832 work_rho, &
833 ener, &
834 v_hartree_gspace=work_v, &
835 calculate_forces=calculate_forces, &
836 itype_of_density="SPIN")
837
838 SELECT CASE (dft_control%sic_method_id)
840 full_scaling = -dft_control%sic_scaling_a
841 CASE (sic_ad)
842 full_scaling = -dft_control%sic_scaling_a*nelec
843 CASE DEFAULT
844 cpabort("Unknown sic method id")
845 END SELECT
846 energy%hartree = energy%hartree + full_scaling*ener
847
848 ! add scaled forces, restoring the old
849 IF (calculate_forces) THEN
850 nforce = 0
851 DO i = 1, SIZE(force)
852 force(i)%ch_pulay(:, :) = force(i)%ch_pulay(:, :)*full_scaling + &
853 store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2))
854 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
855 END DO
856 END IF
857
858 IF (.NOT. just_energy) THEN
859 ALLOCATE (v_sic_rspace)
860 CALL auxbas_pw_pool%create_pw(v_sic_rspace)
861 CALL pw_transfer(work_v, v_sic_rspace)
862 ! also take into account the scaling (in addition to the volume element)
863 CALL pw_scale(v_sic_rspace, &
864 dft_control%sic_scaling_a*v_sic_rspace%pw_grid%dvol)
865 END IF
866
867 CALL auxbas_pw_pool%give_back_pw(work_rho)
868 CALL auxbas_pw_pool%give_back_pw(work_v)
869
870 END SUBROUTINE calc_v_sic_rspace
871
872! **************************************************************************************************
873!> \brief ...
874!> \param qs_env ...
875!> \param rho ...
876! **************************************************************************************************
877 SUBROUTINE print_densities(qs_env, rho)
878 TYPE(qs_environment_type), POINTER :: qs_env
879 TYPE(qs_rho_type), POINTER :: rho
880
881 INTEGER :: img, ispin, n_electrons, output_unit
882 REAL(dp) :: tot1_h, tot1_s, tot_rho_r, trace, &
883 trace_tmp
884 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r_arr
885 TYPE(cell_type), POINTER :: cell
886 TYPE(cp_logger_type), POINTER :: logger
887 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, rho_ao
888 TYPE(dft_control_type), POINTER :: dft_control
889 TYPE(qs_charges_type), POINTER :: qs_charges
890 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
891 TYPE(section_vals_type), POINTER :: input, scf_section
892
893 NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, &
894 dft_control, tot_rho_r_arr, rho_ao)
895
896 logger => cp_get_default_logger()
897
898 CALL get_qs_env(qs_env, &
899 qs_kind_set=qs_kind_set, &
900 cell=cell, qs_charges=qs_charges, &
901 input=input, &
902 matrix_s_kp=matrix_s, &
903 dft_control=dft_control)
904
905 CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
906
907 scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
908 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", &
909 extension=".scfLog")
910
911 CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao)
912 n_electrons = n_electrons - dft_control%charge
913 tot_rho_r = accurate_sum(tot_rho_r_arr)
914
915 trace = 0
916 IF (btest(cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES"), cp_p_file)) THEN
917 DO ispin = 1, dft_control%nspins
918 DO img = 1, dft_control%nimages
919 CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp)
920 trace = trace + trace_tmp
921 END DO
922 END DO
923 END IF
924
925 IF (output_unit > 0) THEN
926 WRITE (unit=output_unit, fmt="(/,T3,A,T41,F20.10)") "Trace(PS):", trace
927 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
928 "Electronic density on regular grids: ", &
929 tot_rho_r, &
930 tot_rho_r + &
931 REAL(n_electrons, dp), &
932 "Core density on regular grids:", &
933 qs_charges%total_rho_core_rspace, &
934 qs_charges%total_rho_core_rspace + &
935 qs_charges%total_rho1_hard_nuc - &
936 REAL(n_electrons + dft_control%charge, dp)
937 END IF
938 IF (dft_control%qs_control%gapw) THEN
939 tot1_h = qs_charges%total_rho1_hard(1)
940 tot1_s = qs_charges%total_rho1_soft(1)
941 DO ispin = 2, dft_control%nspins
942 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
943 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
944 END DO
945 IF (output_unit > 0) THEN
946 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
947 "Hard and soft densities (Lebedev):", &
948 tot1_h, tot1_s
949 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
950 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
951 tot_rho_r + tot1_h - tot1_s, &
952 "Total charge density (r-space): ", &
953 tot_rho_r + tot1_h - tot1_s &
954 + qs_charges%total_rho_core_rspace &
955 + qs_charges%total_rho1_hard_nuc
956 IF (qs_charges%total_rho1_hard_nuc /= 0.0_dp) THEN
957 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
958 "Total CNEO nuc. char. den. (Lebedev): ", &
959 qs_charges%total_rho1_hard_nuc, &
960 "Total CNEO soft char. den. (Lebedev): ", &
961 qs_charges%total_rho1_soft_nuc_lebedev, &
962 "Total CNEO soft char. den. (r-space): ", &
963 qs_charges%total_rho1_soft_nuc_rspace, &
964 "Total soft Rho_e+n+0 (g-space):", &
965 qs_charges%total_rho_gspace
966 ELSE
967 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
968 "Total Rho_soft + Rho0_soft (g-space):", &
969 qs_charges%total_rho_gspace
970 END IF
971 END IF
972 qs_charges%background = tot_rho_r + tot1_h - tot1_s + &
973 qs_charges%total_rho_core_rspace + &
974 qs_charges%total_rho1_hard_nuc
975 ! only add total_rho1_hard_nuc for gapw as cneo requires gapw
976 ELSE IF (dft_control%qs_control%gapw_xc) THEN
977 tot1_h = qs_charges%total_rho1_hard(1)
978 tot1_s = qs_charges%total_rho1_soft(1)
979 DO ispin = 2, dft_control%nspins
980 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
981 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
982 END DO
983 IF (output_unit > 0) THEN
984 WRITE (unit=output_unit, fmt="(/,(T3,A,T41,2F20.10))") &
985 "Hard and soft densities (Lebedev):", &
986 tot1_h, tot1_s
987 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
988 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
989 accurate_sum(tot_rho_r_arr) + tot1_h - tot1_s
990 END IF
991 qs_charges%background = tot_rho_r + &
992 qs_charges%total_rho_core_rspace
993 ELSE
994 IF (output_unit > 0) THEN
995 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
996 "Total charge density on r-space grids: ", &
997 tot_rho_r + &
998 qs_charges%total_rho_core_rspace, &
999 "Total charge density g-space grids: ", &
1000 qs_charges%total_rho_gspace
1001 END IF
1002 qs_charges%background = tot_rho_r + &
1003 qs_charges%total_rho_core_rspace
1004 END IF
1005 IF (output_unit > 0) WRITE (unit=output_unit, fmt="()")
1006 qs_charges%background = qs_charges%background/cell%deth
1007
1008 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
1009 "PRINT%TOTAL_DENSITIES")
1010
1011 END SUBROUTINE print_densities
1012
1013! **************************************************************************************************
1014!> \brief Print detailed energies
1015!>
1016!> \param qs_env ...
1017!> \param dft_control ...
1018!> \param input ...
1019!> \param energy ...
1020!> \param mulliken_order_p ...
1021!> \par History
1022!> refactoring 04.03.2011 [MI]
1023!> \author
1024! **************************************************************************************************
1025 SUBROUTINE print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
1026
1027 TYPE(qs_environment_type), POINTER :: qs_env
1028 TYPE(dft_control_type), POINTER :: dft_control
1029 TYPE(section_vals_type), POINTER :: input
1030 TYPE(qs_energy_type), POINTER :: energy
1031 REAL(kind=dp), INTENT(IN) :: mulliken_order_p
1032
1033 INTEGER :: bc, n, output_unit, psolver
1034 REAL(kind=dp) :: ddapc_order_p, implicit_ps_ehartree, &
1035 s2_order_p
1036 TYPE(cp_logger_type), POINTER :: logger
1037 TYPE(pw_env_type), POINTER :: pw_env
1038
1039 logger => cp_get_default_logger()
1040
1041 NULLIFY (pw_env)
1042 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1043 psolver = pw_env%poisson_env%parameters%solver
1044
1045 output_unit = cp_print_key_unit_nr(logger, input, "DFT%SCF%PRINT%DETAILED_ENERGY", &
1046 extension=".scfLog")
1047 IF (output_unit > 0) THEN
1048 IF (dft_control%do_admm) THEN
1049 WRITE (unit=output_unit, fmt="((T3,A,T60,F20.10))") &
1050 "Wfn fit exchange-correlation energy: ", energy%exc_aux_fit
1051 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1052 WRITE (unit=output_unit, fmt="((T3,A,T60,F20.10))") &
1053 "Wfn fit soft/hard atomic rho1 Exc contribution: ", energy%exc1_aux_fit
1054 END IF
1055 END IF
1056 IF (dft_control%do_admm) THEN
1057 IF (psolver .EQ. pw_poisson_implicit) THEN
1058 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1059 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1060 SELECT CASE (bc)
1062 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1063 "Core Hamiltonian energy: ", energy%core, &
1064 "Hartree energy: ", implicit_ps_ehartree, &
1065 "Electric enthalpy: ", energy%hartree, &
1066 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1067 CASE (periodic_bc, neumann_bc)
1068 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1069 "Core Hamiltonian energy: ", energy%core, &
1070 "Hartree energy: ", energy%hartree, &
1071 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1072 END SELECT
1073 ELSE
1074 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1075 "Core Hamiltonian energy: ", energy%core, &
1076 "Hartree energy: ", energy%hartree, &
1077 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1078 END IF
1079 ELSE
1080!ZMP to print some variables at each step
1081 IF (dft_control%apply_external_density) THEN
1082 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1083 "DOING ZMP CALCULATION FROM EXTERNAL DENSITY "
1084 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1085 "Core Hamiltonian energy: ", energy%core, &
1086 "Hartree energy: ", energy%hartree
1087 ELSE IF (dft_control%apply_external_vxc) THEN
1088 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1089 "DOING ZMP READING EXTERNAL VXC "
1090 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1091 "Core Hamiltonian energy: ", energy%core, &
1092 "Hartree energy: ", energy%hartree
1093 ELSE
1094 IF (psolver .EQ. pw_poisson_implicit) THEN
1095 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1096 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1097 SELECT CASE (bc)
1099 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1100 "Core Hamiltonian energy: ", energy%core, &
1101 "Hartree energy: ", implicit_ps_ehartree, &
1102 "Electric enthalpy: ", energy%hartree, &
1103 "Exchange-correlation energy: ", energy%exc
1104 CASE (periodic_bc, neumann_bc)
1105 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1106 "Core Hamiltonian energy: ", energy%core, &
1107 "Hartree energy: ", energy%hartree, &
1108 "Exchange-correlation energy: ", energy%exc
1109 END SELECT
1110 ELSE
1111 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1112 "Core Hamiltonian energy: ", energy%core, &
1113 "Hartree energy: ", energy%hartree, &
1114 "Exchange-correlation energy: ", energy%exc
1115 END IF
1116 END IF
1117 END IF
1118
1119 IF (dft_control%apply_external_density) THEN
1120 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1121 "Integral of the (density * v_xc): ", energy%exc
1122 END IF
1123
1124 IF (energy%e_hartree /= 0.0_dp) &
1125 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1126 "Coulomb (electron-electron) energy: ", energy%e_hartree
1127 IF (energy%dispersion /= 0.0_dp) &
1128 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1129 "Dispersion energy: ", energy%dispersion
1130 IF (energy%efield /= 0.0_dp) &
1131 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1132 "Electric field interaction energy: ", energy%efield
1133 IF (energy%gcp /= 0.0_dp) &
1134 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1135 "gCP energy: ", energy%gcp
1136 IF (dft_control%qs_control%gapw) THEN
1137 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1138 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit, &
1139 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
1140 END IF
1141 IF (dft_control%qs_control%gapw_xc) THEN
1142 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1143 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit
1144 END IF
1145 IF (dft_control%dft_plus_u) THEN
1146 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1147 "DFT+U energy:", energy%dft_plus_u
1148 END IF
1149 IF (qs_env%qmmm) THEN
1150 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1151 "QM/MM Electrostatic energy: ", energy%qmmm_el
1152 IF (qs_env%qmmm_env_qm%image_charge) THEN
1153 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1154 "QM/MM image charge energy: ", energy%image_charge
1155 END IF
1156 END IF
1157 IF (dft_control%qs_control%mulliken_restraint) THEN
1158 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1159 "Mulliken restraint (order_p,energy) : ", mulliken_order_p, energy%mulliken
1160 END IF
1161 IF (dft_control%qs_control%ddapc_restraint) THEN
1162 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
1163 ddapc_order_p = &
1164 dft_control%qs_control%ddapc_restraint_control(n)%ddapc_order_p
1165 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1166 "DDAPC restraint (order_p,energy) : ", ddapc_order_p, energy%ddapc_restraint(n)
1167 END DO
1168 END IF
1169 IF (dft_control%qs_control%s2_restraint) THEN
1170 s2_order_p = dft_control%qs_control%s2_restraint_control%s2_order_p
1171 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1172 "S2 restraint (order_p,energy) : ", s2_order_p, energy%s2_restraint
1173 END IF
1174 IF (energy%core_cneo /= 0.0_dp) THEN
1175 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1176 "CNEO| quantum nuclear core energy: ", energy%core_cneo
1177 END IF
1178
1179 END IF ! output_unit
1180 CALL cp_print_key_finished_output(output_unit, logger, input, &
1181 "DFT%SCF%PRINT%DETAILED_ENERGY")
1182
1183 END SUBROUTINE print_detailed_energy
1184
1185! **************************************************************************************************
1186!> \brief compute matrix_vxc, defined via the potential created by qs_vxc_create
1187!> ignores things like tau functional, gapw, sic, ...
1188!> so only OK for GGA & GPW right now
1189!> \param qs_env ...
1190!> \param v_rspace ...
1191!> \param matrix_vxc ...
1192!> \par History
1193!> created 23.10.2012 [Joost VandeVondele]
1194!> \author
1195! **************************************************************************************************
1196 SUBROUTINE compute_matrix_vxc(qs_env, v_rspace, matrix_vxc)
1197 TYPE(qs_environment_type), POINTER :: qs_env
1198 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
1199 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
1200
1201 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_matrix_vxc'
1202
1203 INTEGER :: handle, ispin
1204 LOGICAL :: gapw
1205 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1206 TYPE(dft_control_type), POINTER :: dft_control
1207
1208 CALL timeset(routinen, handle)
1209
1210 ! create the matrix using matrix_ks as a template
1211 IF (ASSOCIATED(matrix_vxc)) THEN
1212 CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1213 END IF
1214 CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
1215 ALLOCATE (matrix_vxc(SIZE(matrix_ks)))
1216 DO ispin = 1, SIZE(matrix_ks)
1217 NULLIFY (matrix_vxc(ispin)%matrix)
1218 CALL dbcsr_init_p(matrix_vxc(ispin)%matrix)
1219 CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, &
1220 name="Matrix VXC of spin "//cp_to_string(ispin))
1221 CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp)
1222 END DO
1223
1224 ! and integrate
1225 CALL get_qs_env(qs_env, dft_control=dft_control)
1226 gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1227 DO ispin = 1, SIZE(matrix_ks)
1228 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1229 hmat=matrix_vxc(ispin), &
1230 qs_env=qs_env, &
1231 calculate_forces=.false., &
1232 gapw=gapw)
1233 ! scale by the volume element... should really become part of integrate_v_rspace
1234 CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw_grid%dvol)
1235 END DO
1236
1237 CALL timestop(handle)
1238
1239 END SUBROUTINE compute_matrix_vxc
1240
1241! **************************************************************************************************
1242!> \brief Sum up all potentials defined on the grid and integrate
1243!>
1244!> \param qs_env ...
1245!> \param ks_matrix ...
1246!> \param rho ...
1247!> \param my_rho ...
1248!> \param vppl_rspace ...
1249!> \param v_rspace_new ...
1250!> \param v_rspace_new_aux_fit ...
1251!> \param v_tau_rspace ...
1252!> \param v_tau_rspace_aux_fit ...
1253!> \param v_sic_rspace ...
1254!> \param v_spin_ddapc_rest_r ...
1255!> \param v_sccs_rspace ...
1256!> \param v_rspace_embed ...
1257!> \param cdft_control ...
1258!> \param calculate_forces ...
1259!> \par History
1260!> - refactoring 04.03.2011 [MI]
1261!> - SCCS implementation (16.10.2013,MK)
1262!> \author
1263! **************************************************************************************************
1264 SUBROUTINE sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, &
1265 vppl_rspace, v_rspace_new, &
1266 v_rspace_new_aux_fit, v_tau_rspace, &
1267 v_tau_rspace_aux_fit, &
1268 v_sic_rspace, v_spin_ddapc_rest_r, &
1269 v_sccs_rspace, v_rspace_embed, cdft_control, &
1270 calculate_forces)
1271
1272 TYPE(qs_environment_type), POINTER :: qs_env
1273 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
1274 TYPE(qs_rho_type), POINTER :: rho
1275 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: my_rho
1276 TYPE(pw_r3d_rs_type), POINTER :: vppl_rspace
1277 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_rspace_new_aux_fit, &
1278 v_tau_rspace, v_tau_rspace_aux_fit
1279 TYPE(pw_r3d_rs_type), POINTER :: v_sic_rspace, v_spin_ddapc_rest_r, &
1280 v_sccs_rspace
1281 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_embed
1282 TYPE(cdft_control_type), POINTER :: cdft_control
1283 LOGICAL, INTENT(in) :: calculate_forces
1284
1285 CHARACTER(LEN=*), PARAMETER :: routinen = 'sum_up_and_integrate'
1286
1287 CHARACTER(LEN=default_string_length) :: basis_type
1288 INTEGER :: handle, igroup, ikind, img, ispin, &
1289 nkind, nspins
1290 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1291 LOGICAL :: do_ppl, gapw, gapw_xc, lrigpw, rigpw
1292 REAL(kind=dp) :: csign, dvol, fadm
1293 TYPE(admm_type), POINTER :: admm_env
1294 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1295 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao, rho_ao_nokp, smat
1296 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, &
1297 matrix_ks_aux_fit_dft, rho_ao_aux, &
1298 rho_ao_kp
1299 TYPE(dft_control_type), POINTER :: dft_control
1300 TYPE(kpoint_type), POINTER :: kpoints
1301 TYPE(lri_density_type), POINTER :: lri_density
1302 TYPE(lri_environment_type), POINTER :: lri_env
1303 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
1304 TYPE(mp_para_env_type), POINTER :: para_env
1305 TYPE(pw_env_type), POINTER :: pw_env
1306 TYPE(pw_poisson_type), POINTER :: poisson_env
1307 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1308 TYPE(pw_r3d_rs_type), POINTER :: v_rspace, vee
1309 TYPE(qs_ks_env_type), POINTER :: ks_env
1310 TYPE(qs_rho_type), POINTER :: rho_aux_fit
1311 TYPE(task_list_type), POINTER :: task_list
1312
1313 CALL timeset(routinen, handle)
1314
1315 NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, &
1316 v_rspace, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, &
1317 ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, &
1318 rho_ao_nokp, ks_env, admm_env, task_list)
1319
1320 CALL get_qs_env(qs_env, &
1321 dft_control=dft_control, &
1322 pw_env=pw_env, &
1323 v_hartree_rspace=v_rspace, &
1324 vee=vee)
1325
1326 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1327 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
1328 gapw = dft_control%qs_control%gapw
1329 gapw_xc = dft_control%qs_control%gapw_xc
1330 do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
1331
1332 rigpw = dft_control%qs_control%rigpw
1333 lrigpw = dft_control%qs_control%lrigpw
1334 IF (lrigpw .OR. rigpw) THEN
1335 CALL get_qs_env(qs_env, &
1336 lri_env=lri_env, &
1337 lri_density=lri_density, &
1338 atomic_kind_set=atomic_kind_set)
1339 END IF
1340
1341 nspins = dft_control%nspins
1342
1343 ! sum up potentials and integrate
1344 IF (ASSOCIATED(v_rspace_new)) THEN
1345 DO ispin = 1, nspins
1346 IF (gapw_xc) THEN
1347 ! SIC not implemented (or at least not tested)
1348 cpassert(dft_control%sic_method_id == sic_none)
1349 !Only the xc potential, because it has to be integrated with the soft basis
1350 CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
1351
1352 ! add the xc part due to v_rspace soft
1353 rho_ao => rho_ao_kp(ispin, :)
1354 ksmat => ks_matrix(ispin, :)
1355 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1356 pmat_kp=rho_ao, hmat_kp=ksmat, &
1357 qs_env=qs_env, &
1358 calculate_forces=calculate_forces, &
1359 gapw=gapw_xc)
1360
1361 ! Now the Hartree potential to be integrated with the full basis
1362 CALL pw_copy(v_rspace, v_rspace_new(ispin))
1363 ELSE
1364 ! Add v_hartree + v_xc = v_rspace_new
1365 CALL pw_axpy(v_rspace, v_rspace_new(ispin), 1.0_dp, v_rspace_new(ispin)%pw_grid%dvol)
1366 END IF ! gapw_xc
1367 IF (dft_control%qs_control%ddapc_explicit_potential) THEN
1368 IF (dft_control%qs_control%ddapc_restraint_is_spin) THEN
1369 IF (ispin == 1) THEN
1370 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1371 ELSE
1372 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), -1.0_dp)
1373 END IF
1374 ELSE
1375 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1376 END IF
1377 END IF
1378 ! CDFT constraint contribution
1379 IF (dft_control%qs_control%cdft) THEN
1380 DO igroup = 1, SIZE(cdft_control%group)
1381 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1383 csign = 1.0_dp
1385 IF (ispin == 1) THEN
1386 csign = 1.0_dp
1387 ELSE
1388 csign = -1.0_dp
1389 END IF
1391 csign = 1.0_dp
1392 IF (ispin == 2) cycle
1394 csign = 1.0_dp
1395 IF (ispin == 1) cycle
1396 CASE DEFAULT
1397 cpabort("Unknown constraint type.")
1398 END SELECT
1399 CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_new(ispin), &
1400 csign*cdft_control%strength(igroup))
1401 END DO
1402 END IF
1403 ! functional derivative of the Hartree energy wrt the density in the presence of dielectric
1404 ! (vhartree + v_eps); v_eps is nonzero only if the dielectric constant is defind as a function
1405 ! of the charge density
1406 IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
1407 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1408 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_new(ispin), dvol)
1409 END IF
1410 ! Add SCCS contribution
1411 IF (dft_control%do_sccs) THEN
1412 CALL pw_axpy(v_sccs_rspace, v_rspace_new(ispin))
1413 END IF
1414 ! External electrostatic potential
1415 IF (dft_control%apply_external_potential) THEN
1416 CALL qmmm_modify_hartree_pot(v_hartree=v_rspace_new(ispin), &
1417 v_qmmm=vee, scale=-1.0_dp)
1418 END IF
1419 IF (do_ppl) THEN
1420 cpassert(.NOT. gapw)
1421 CALL pw_axpy(vppl_rspace, v_rspace_new(ispin), vppl_rspace%pw_grid%dvol)
1422 END IF
1423 ! the electrostatic sic contribution
1424 SELECT CASE (dft_control%sic_method_id)
1425 CASE (sic_none)
1426 !
1428 IF (ispin == 1) THEN
1429 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1430 ELSE
1431 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), 1.0_dp)
1432 END IF
1433 CASE (sic_ad)
1434 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1435 CASE (sic_eo)
1436 ! NOTHING TO BE DONE
1437 END SELECT
1438 ! DFT embedding
1439 IF (dft_control%apply_embed_pot) THEN
1440 CALL pw_axpy(v_rspace_embed(ispin), v_rspace_new(ispin), v_rspace_embed(ispin)%pw_grid%dvol)
1441 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1442 END IF
1443 IF (lrigpw) THEN
1444 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1445 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1446 DO ikind = 1, nkind
1447 lri_v_int(ikind)%v_int = 0.0_dp
1448 END DO
1449 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1450 lri_v_int, calculate_forces, "LRI_AUX")
1451 DO ikind = 1, nkind
1452 CALL para_env%sum(lri_v_int(ikind)%v_int)
1453 END DO
1454 IF (lri_env%exact_1c_terms) THEN
1455 rho_ao => my_rho(ispin, :)
1456 ksmat => ks_matrix(ispin, :)
1457 CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), ksmat(1)%matrix, &
1458 rho_ao(1)%matrix, qs_env, &
1459 calculate_forces, "ORB")
1460 END IF
1461 IF (lri_env%ppl_ri) THEN
1462 CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1463 END IF
1464 ELSEIF (rigpw) THEN
1465 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1466 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1467 DO ikind = 1, nkind
1468 lri_v_int(ikind)%v_int = 0.0_dp
1469 END DO
1470 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1471 lri_v_int, calculate_forces, "RI_HXC")
1472 DO ikind = 1, nkind
1473 CALL para_env%sum(lri_v_int(ikind)%v_int)
1474 END DO
1475 ELSE
1476 rho_ao => my_rho(ispin, :)
1477 ksmat => ks_matrix(ispin, :)
1478 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1479 pmat_kp=rho_ao, hmat_kp=ksmat, &
1480 qs_env=qs_env, &
1481 calculate_forces=calculate_forces, &
1482 gapw=gapw)
1483 END IF
1484 CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
1485 END DO ! ispin
1486
1487 SELECT CASE (dft_control%sic_method_id)
1488 CASE (sic_none)
1490 CALL auxbas_pw_pool%give_back_pw(v_sic_rspace)
1491 DEALLOCATE (v_sic_rspace)
1492 END SELECT
1493 DEALLOCATE (v_rspace_new)
1494
1495 ELSE
1496 ! not implemented (or at least not tested)
1497 cpassert(dft_control%sic_method_id == sic_none)
1498 cpassert(.NOT. dft_control%qs_control%ddapc_restraint_is_spin)
1499 DO ispin = 1, nspins
1500 ! extra contribution attributed to the dependency of the dielectric constant to the charge density
1501 IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
1502 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1503 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace, dvol)
1504 END IF
1505 ! Add SCCS contribution
1506 IF (dft_control%do_sccs) THEN
1507 CALL pw_axpy(v_sccs_rspace, v_rspace)
1508 END IF
1509 ! DFT embedding
1510 IF (dft_control%apply_embed_pot) THEN
1511 CALL pw_axpy(v_rspace_embed(ispin), v_rspace, v_rspace_embed(ispin)%pw_grid%dvol)
1512 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1513 END IF
1514 IF (lrigpw) THEN
1515 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1516 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1517 DO ikind = 1, nkind
1518 lri_v_int(ikind)%v_int = 0.0_dp
1519 END DO
1520 CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
1521 lri_v_int, calculate_forces, "LRI_AUX")
1522 DO ikind = 1, nkind
1523 CALL para_env%sum(lri_v_int(ikind)%v_int)
1524 END DO
1525 IF (lri_env%exact_1c_terms) THEN
1526 rho_ao => my_rho(ispin, :)
1527 ksmat => ks_matrix(ispin, :)
1528 CALL integrate_v_rspace_diagonal(v_rspace, ksmat(1)%matrix, &
1529 rho_ao(1)%matrix, qs_env, &
1530 calculate_forces, "ORB")
1531 END IF
1532 IF (lri_env%ppl_ri) THEN
1533 CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1534 END IF
1535 ELSEIF (rigpw) THEN
1536 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1537 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1538 DO ikind = 1, nkind
1539 lri_v_int(ikind)%v_int = 0.0_dp
1540 END DO
1541 CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
1542 lri_v_int, calculate_forces, "RI_HXC")
1543 DO ikind = 1, nkind
1544 CALL para_env%sum(lri_v_int(ikind)%v_int)
1545 END DO
1546 ELSE
1547 rho_ao => my_rho(ispin, :)
1548 ksmat => ks_matrix(ispin, :)
1549 CALL integrate_v_rspace(v_rspace=v_rspace, &
1550 pmat_kp=rho_ao, &
1551 hmat_kp=ksmat, &
1552 qs_env=qs_env, &
1553 calculate_forces=calculate_forces, &
1554 gapw=gapw)
1555 END IF
1556 END DO
1557 END IF ! ASSOCIATED(v_rspace_new)
1558
1559 ! **** LRIGPW: KS matrix from integrated potential
1560 IF (lrigpw) THEN
1561 CALL get_qs_env(qs_env, ks_env=ks_env)
1562 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1563 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1564 DO ispin = 1, nspins
1565 ksmat => ks_matrix(ispin, :)
1566 CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set, &
1567 cell_to_index=cell_to_index)
1568 END DO
1569 IF (calculate_forces) THEN
1570 CALL calculate_lri_forces(lri_env, lri_density, qs_env, rho_ao_kp, atomic_kind_set)
1571 END IF
1572 ELSEIF (rigpw) THEN
1573 CALL get_qs_env(qs_env, matrix_s=smat)
1574 DO ispin = 1, nspins
1575 CALL calculate_ri_ks_matrix(lri_env, lri_v_int, ks_matrix(ispin, 1)%matrix, &
1576 smat(1)%matrix, atomic_kind_set, ispin)
1577 END DO
1578 IF (calculate_forces) THEN
1579 rho_ao_nokp => rho_ao_kp(:, 1)
1580 CALL calculate_ri_forces(lri_env, lri_density, qs_env, rho_ao_nokp, atomic_kind_set)
1581 END IF
1582 END IF
1583
1584 IF (ASSOCIATED(v_tau_rspace)) THEN
1585 IF (lrigpw .OR. rigpw) THEN
1586 cpabort("LRIGPW/RIGPW not implemented for meta-GGAs")
1587 END IF
1588 DO ispin = 1, nspins
1589 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1590
1591 rho_ao => rho_ao_kp(ispin, :)
1592 ksmat => ks_matrix(ispin, :)
1593 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1594 pmat_kp=rho_ao, hmat_kp=ksmat, &
1595 qs_env=qs_env, &
1596 calculate_forces=calculate_forces, compute_tau=.true., &
1597 gapw=gapw)
1598 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1599 END DO
1600 DEALLOCATE (v_tau_rspace)
1601 END IF
1602
1603 ! Add contributions from ADMM if requested
1604 IF (dft_control%do_admm) THEN
1605 CALL get_qs_env(qs_env, admm_env=admm_env)
1606 CALL get_admm_env(admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
1607 matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft)
1608 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
1609 IF (ASSOCIATED(v_rspace_new_aux_fit)) THEN
1610 DO ispin = 1, nspins
1611 ! Calculate the xc potential
1612 CALL pw_scale(v_rspace_new_aux_fit(ispin), v_rspace_new_aux_fit(ispin)%pw_grid%dvol)
1613
1614 ! set matrix_ks_aux_fit_dft = matrix_ks_aux_fit(k_HF)
1615 DO img = 1, dft_control%nimages
1616 CALL dbcsr_copy(matrix_ks_aux_fit_dft(ispin, img)%matrix, matrix_ks_aux_fit(ispin, img)%matrix, &
1617 name="DFT exch. part of matrix_ks_aux_fit")
1618 END DO
1619
1620 ! Add potential to ks_matrix aux_fit, skip integration if no DFT correction
1621
1622 IF (admm_env%aux_exch_func .NE. do_admm_aux_exch_func_none) THEN
1623
1624 !GPW by default. IF GAPW, then take relevant task list and basis
1625 CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
1626 basis_type = "AUX_FIT"
1627 IF (admm_env%do_gapw) THEN
1628 task_list => admm_env%admm_gapw_env%task_list
1629 basis_type = "AUX_FIT_SOFT"
1630 END IF
1631 fadm = 1.0_dp
1632 ! Calculate bare scaling of force according to Merlot, 1. IF: ADMMP, 2. IF: ADMMS,
1633 IF (admm_env%do_admmp) THEN
1634 fadm = admm_env%gsi(ispin)**2
1635 ELSE IF (admm_env%do_admms) THEN
1636 fadm = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)
1637 END IF
1638
1639 rho_ao => rho_ao_aux(ispin, :)
1640 ksmat => matrix_ks_aux_fit(ispin, :)
1641
1642 CALL integrate_v_rspace(v_rspace=v_rspace_new_aux_fit(ispin), &
1643 pmat_kp=rho_ao, &
1644 hmat_kp=ksmat, &
1645 qs_env=qs_env, &
1646 calculate_forces=calculate_forces, &
1647 force_adm=fadm, &
1648 gapw=.false., & !even if actual GAPW calculation, want to use AUX_FIT_SOFT
1649 basis_type=basis_type, &
1650 task_list_external=task_list)
1651 END IF
1652
1653 ! matrix_ks_aux_fit_dft(x_DFT)=matrix_ks_aux_fit_dft(old,k_HF)-matrix_ks_aux_fit(k_HF-x_DFT)
1654 DO img = 1, dft_control%nimages
1655 CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, &
1656 matrix_ks_aux_fit(ispin, img)%matrix, 1.0_dp, -1.0_dp)
1657 END DO
1658
1659 CALL auxbas_pw_pool%give_back_pw(v_rspace_new_aux_fit(ispin))
1660 END DO
1661 DEALLOCATE (v_rspace_new_aux_fit)
1662 END IF
1663 ! Clean up v_tau_rspace_aux_fit, which is actually not needed
1664 IF (ASSOCIATED(v_tau_rspace_aux_fit)) THEN
1665 DO ispin = 1, nspins
1666 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace_aux_fit(ispin))
1667 END DO
1668 DEALLOCATE (v_tau_rspace_aux_fit)
1669 END IF
1670 END IF
1671
1672 IF (dft_control%apply_embed_pot) DEALLOCATE (v_rspace_embed)
1673
1674 CALL timestop(handle)
1675
1676 END SUBROUTINE sum_up_and_integrate
1677
1678!**************************************************************************
1679!> \brief Calculate the ZMP potential and energy as in Zhao, Morrison Parr
1680!> PRA 50i, 2138 (1994)
1681!> V_c^\lambda defined as int_rho-rho_0/r-r' or rho-rho_0 times a Lagrange
1682!> multiplier, plus Fermi-Amaldi potential that should give the V_xc in the
1683!> limit \lambda --> \infty
1684!>
1685!> \param qs_env ...
1686!> \param v_rspace_new ...
1687!> \param rho ...
1688!> \param exc ...
1689!> \author D. Varsano [daniele.varsano@nano.cnr.it]
1690! **************************************************************************************************
1691 SUBROUTINE calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
1692
1693 TYPE(qs_environment_type), POINTER :: qs_env
1694 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new
1695 TYPE(qs_rho_type), POINTER :: rho
1696 REAL(kind=dp) :: exc
1697
1698 CHARACTER(*), PARAMETER :: routinen = 'calculate_zmp_potential'
1699
1700 INTEGER :: handle, my_val, nelectron, nspins
1701 INTEGER, DIMENSION(2) :: nelectron_spin
1702 LOGICAL :: do_zmp_read, fermi_amaldi
1703 REAL(kind=dp) :: lambda
1704 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_ext_r
1705 TYPE(dft_control_type), POINTER :: dft_control
1706 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ext_g, rho_g
1707 TYPE(pw_env_type), POINTER :: pw_env
1708 TYPE(pw_poisson_type), POINTER :: poisson_env
1709 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1710 TYPE(pw_r3d_rs_type) :: v_xc_rspace
1711 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1712 TYPE(qs_ks_env_type), POINTER :: ks_env
1713 TYPE(section_vals_type), POINTER :: ext_den_section, input
1714
1715!, v_h_gspace, &
1716
1717 CALL timeset(routinen, handle)
1718 NULLIFY (auxbas_pw_pool)
1719 NULLIFY (pw_env)
1720 NULLIFY (poisson_env)
1721 NULLIFY (v_rspace_new)
1722 NULLIFY (dft_control)
1723 NULLIFY (rho_r, rho_g, tot_rho_ext_r, rho_ext_g)
1724 CALL get_qs_env(qs_env=qs_env, &
1725 pw_env=pw_env, &
1726 ks_env=ks_env, &
1727 rho=rho, &
1728 input=input, &
1729 nelectron_spin=nelectron_spin, &
1730 dft_control=dft_control)
1731 CALL pw_env_get(pw_env=pw_env, &
1732 auxbas_pw_pool=auxbas_pw_pool, &
1733 poisson_env=poisson_env)
1734 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
1735 nspins = 1
1736 ALLOCATE (v_rspace_new(nspins))
1737 CALL auxbas_pw_pool%create_pw(pw=v_rspace_new(1))
1738 CALL auxbas_pw_pool%create_pw(pw=v_xc_rspace)
1739
1740 CALL pw_zero(v_rspace_new(1))
1741 do_zmp_read = dft_control%apply_external_vxc
1742 IF (do_zmp_read) THEN
1743 CALL pw_copy(qs_env%external_vxc, v_rspace_new(1))
1744 exc = accurate_dot_product(v_rspace_new(1)%array, rho_r(1)%array)* &
1745 v_rspace_new(1)%pw_grid%dvol
1746 ELSE
1747 block
1748 REAL(kind=dp) :: factor
1749 TYPE(pw_c1d_gs_type) :: rho_eff_gspace, v_xc_gspace
1750 CALL auxbas_pw_pool%create_pw(pw=rho_eff_gspace)
1751 CALL auxbas_pw_pool%create_pw(pw=v_xc_gspace)
1752 CALL pw_zero(rho_eff_gspace)
1753 CALL pw_zero(v_xc_gspace)
1754 CALL pw_zero(v_xc_rspace)
1755 factor = pw_integrate_function(rho_g(1))
1756 CALL qs_rho_get(qs_env%rho_external, &
1757 rho_g=rho_ext_g, &
1758 tot_rho_r=tot_rho_ext_r)
1759 factor = tot_rho_ext_r(1)/factor
1760
1761 CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1762 CALL pw_axpy(rho_ext_g(1), rho_eff_gspace, alpha=-1.0_dp)
1763 ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
1764 CALL section_vals_val_get(ext_den_section, "LAMBDA", r_val=lambda)
1765 CALL section_vals_val_get(ext_den_section, "ZMP_CONSTRAINT", i_val=my_val)
1766 CALL section_vals_val_get(ext_den_section, "FERMI_AMALDI", l_val=fermi_amaldi)
1767
1768 CALL pw_scale(rho_eff_gspace, a=lambda)
1769 nelectron = nelectron_spin(1)
1770 factor = -1.0_dp/nelectron
1771 CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1772
1773 CALL pw_poisson_solve(poisson_env, rho_eff_gspace, vhartree=v_xc_gspace)
1774 CALL pw_transfer(v_xc_gspace, v_rspace_new(1))
1775 CALL pw_copy(v_rspace_new(1), v_xc_rspace)
1776
1777 exc = 0.0_dp
1778 exc = pw_integral_ab(v_rspace_new(1), rho_r(1))
1779
1780!Note that this is not the xc energy but \int(\rho*v_xc)
1781!Vxc---> v_rspace_new
1782!Exc---> energy%exc
1783 CALL auxbas_pw_pool%give_back_pw(rho_eff_gspace)
1784 CALL auxbas_pw_pool%give_back_pw(v_xc_gspace)
1785 END block
1786 END IF
1787
1788 CALL auxbas_pw_pool%give_back_pw(v_xc_rspace)
1789
1790 CALL timestop(handle)
1791
1792 END SUBROUTINE calculate_zmp_potential
1793
1794! **************************************************************************************************
1795!> \brief ...
1796!> \param qs_env ...
1797!> \param rho ...
1798!> \param v_rspace_embed ...
1799!> \param dft_control ...
1800!> \param embed_corr ...
1801!> \param just_energy ...
1802! **************************************************************************************************
1803 SUBROUTINE get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, &
1804 just_energy)
1805 TYPE(qs_environment_type), POINTER :: qs_env
1806 TYPE(qs_rho_type), POINTER :: rho
1807 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_embed
1808 TYPE(dft_control_type), POINTER :: dft_control
1809 REAL(kind=dp) :: embed_corr
1810 LOGICAL :: just_energy
1811
1812 CHARACTER(*), PARAMETER :: routinen = 'get_embed_potential_energy'
1813
1814 INTEGER :: handle, ispin
1815 REAL(kind=dp) :: embed_corr_local
1816 TYPE(pw_env_type), POINTER :: pw_env
1817 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1818 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1819
1820 CALL timeset(routinen, handle)
1821
1822 NULLIFY (auxbas_pw_pool)
1823 NULLIFY (pw_env)
1824 NULLIFY (rho_r)
1825 CALL get_qs_env(qs_env=qs_env, &
1826 pw_env=pw_env, &
1827 rho=rho)
1828 CALL pw_env_get(pw_env=pw_env, &
1829 auxbas_pw_pool=auxbas_pw_pool)
1830 CALL qs_rho_get(rho, rho_r=rho_r)
1831 ALLOCATE (v_rspace_embed(dft_control%nspins))
1832
1833 embed_corr = 0.0_dp
1834
1835 DO ispin = 1, dft_control%nspins
1836 CALL auxbas_pw_pool%create_pw(pw=v_rspace_embed(ispin))
1837 CALL pw_zero(v_rspace_embed(ispin))
1838
1839 CALL pw_copy(qs_env%embed_pot, v_rspace_embed(ispin))
1840 embed_corr_local = 0.0_dp
1841
1842 ! Spin embedding potential in open-shell case
1843 IF (dft_control%nspins .EQ. 2) THEN
1844 IF (ispin .EQ. 1) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), 1.0_dp)
1845 IF (ispin .EQ. 2) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), -1.0_dp)
1846 END IF
1847 ! Integrate the density*potential
1848 embed_corr_local = pw_integral_ab(v_rspace_embed(ispin), rho_r(ispin))
1849
1850 embed_corr = embed_corr + embed_corr_local
1851
1852 END DO
1853
1854 ! If only energy requiested we delete the potential
1855 IF (just_energy) THEN
1856 DO ispin = 1, dft_control%nspins
1857 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1858 END DO
1859 DEALLOCATE (v_rspace_embed)
1860 END IF
1861
1862 CALL timestop(handle)
1863
1864 END SUBROUTINE get_embed_potential_energy
1865
1866END MODULE qs_ks_utils
Types and set/get functions for auxiliary density matrix methods.
Definition admm_types.F:15
subroutine, public get_admm_env(admm_env, mo_derivs_aux_fit, mos_aux_fit, sab_aux_fit, sab_aux_fit_asymm, sab_aux_fit_vs_orb, matrix_s_aux_fit, matrix_s_aux_fit_kp, matrix_s_aux_fit_vs_orb, matrix_s_aux_fit_vs_orb_kp, task_list_aux_fit, matrix_ks_aux_fit, matrix_ks_aux_fit_kp, matrix_ks_aux_fit_im, matrix_ks_aux_fit_dft, matrix_ks_aux_fit_hfx, matrix_ks_aux_fit_dft_kp, matrix_ks_aux_fit_hfx_kp, rho_aux_fit, rho_aux_fit_buffer, admm_dm)
Get routine for the ADMM env.
Definition admm_types.F:593
Define the atomic kind types and their sub types.
Handles all functions related to the CELL.
Definition cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_release_p(matrix)
...
subroutine, public dbcsr_scale(matrix, alpha_scalar)
...
subroutine, public dbcsr_deallocate_matrix(matrix)
...
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
subroutine, public dbcsr_multiply(transa, transb, alpha, matrix_a, matrix_b, beta, matrix_c, first_row, last_row, first_column, last_column, first_k, last_k, retain_sparsity, filter_eps, flop)
...
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
subroutine, public dbcsr_init_p(matrix)
...
subroutine, public dbcsr_set(matrix, alpha)
...
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
subroutine, public dbcsr_dot(matrix_a, matrix_b, trace)
Computes the dot product of two matrices, also known as the trace of their matrix product.
subroutine, public dbcsr_scale_by_vector(matrix, alpha, side)
Scales the rows/columns of given matrix.
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
subroutine, public copy_fm_to_dbcsr(fm, matrix, keep_sparsity)
Copy a BLACS matrix to a dbcsr matrix.
Density Derived atomic point charges from a QM calculation (see Bloechl, J. Chem. Phys....
Definition cp_ddapc.F:15
subroutine, public cp_ddapc_apply_cd(qs_env, rho_tot_gspace, energy, v_hartree_gspace, calculate_forces, itype_of_density)
Routine to couple/decouple periodic images with the Bloechl scheme.
Definition cp_ddapc.F:226
represent the structure of a full matrix
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Utilities for hfx and admm methods.
subroutine, public tddft_hfx_matrix(matrix_ks, rho_ao, qs_env, update_energy, recalc_integrals, external_hfx_sections, external_x_data, external_para_env)
Add the hfx contributions to the Hamiltonian.
Routines to calculate derivatives with respect to basis function origin.
subroutine, public derivatives_four_center(qs_env, rho_ao, rho_ao_resp, hfx_section, para_env, irep, use_virial, adiabatic_rescale_factor, resp_only, external_x_data)
computes four center derivatives for a full basis set and updates the forcesfock_4c arrays....
Types and set/get functions for HFX.
Definition hfx_types.F:15
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public sic_list_unpaired
integer, parameter, public sic_mauri_spz
integer, parameter, public cdft_beta_constraint
integer, parameter, public cdft_magnetization_constraint
integer, parameter, public sic_list_all
integer, parameter, public cdft_charge_constraint
integer, parameter, public sic_eo
integer, parameter, public do_ppl_grid
integer, parameter, public do_admm_aux_exch_func_none
integer, parameter, public sic_mauri_us
integer, parameter, public sic_none
integer, parameter, public cdft_alpha_constraint
integer, parameter, public sic_ad
objects that represent the structure of input sections and the data contained in an input section
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
integer, parameter, public default_string_length
Definition kinds.F:57
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Calculates integral matrices for LRIGPW method lri : local resolution of the identity.
subroutine, public v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
...
contains the types and subroutines for dealing with the lri_env lri : local resolution of the identit...
Calculates forces for LRIGPW method lri : local resolution of the identity.
Definition lri_forces.F:16
subroutine, public calculate_lri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the lri forces
Definition lri_forces.F:81
subroutine, public calculate_ri_forces(lri_env, lri_density, qs_env, pmatrix, atomic_kind_set)
calculates the ri forces
Definition lri_forces.F:673
routines that build the Kohn-Sham matrix for the LRIGPW and xc parts
subroutine, public calculate_lri_ks_matrix(lri_env, lri_v_int, h_matrix, atomic_kind_set, cell_to_index)
update of LRIGPW KS matrix
subroutine, public calculate_ri_ks_matrix(lri_env, lri_v_int, h_matrix, s_matrix, atomic_kind_set, ispin)
update of RIGPW KS matrix
Interface to the message passing library MPI.
Types containing essential information for running implicit (iterative) Poisson solver.
integer, parameter, public neumann_bc
integer, parameter, public mixed_bc
integer, parameter, public mixed_periodic_bc
integer, parameter, public periodic_bc
container for various plainwaves related things
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
functions related to the poisson solver on regular grids
integer, parameter, public pw_poisson_implicit
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Defines CDFT control structures.
container for information about total charges on the grids
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr, npgf_seg, cneo_potential_present, nkind_q, natom_q)
Get attributes of an atomic kind set.
subroutine, public qmmm_modify_hartree_pot(v_hartree, v_qmmm, scale)
Modify the hartree potential in order to include the QM/MM correction.
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_utils.F:22
subroutine, public print_densities(qs_env, rho)
...
subroutine, public get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, just_energy)
...
subroutine, public low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, calculate_forces, auxbas_pw_pool)
do ROKS calculations yielding low spin states
subroutine, public sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, vppl_rspace, v_rspace_new, v_rspace_new_aux_fit, v_tau_rspace, v_tau_rspace_aux_fit, v_sic_rspace, v_spin_ddapc_rest_r, v_sccs_rspace, v_rspace_embed, cdft_control, calculate_forces)
Sum up all potentials defined on the grid and integrate.
subroutine, public print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
Print detailed energies.
subroutine, public calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
Calculate the ZMP potential and energy as in Zhao, Morrison Parr PRA 50i, 2138 (1994) V_c^\lambda def...
subroutine, public calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, just_energy, calculate_forces, auxbas_pw_pool)
do sic calculations on the spin density
subroutine, public sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, calculate_forces, auxbas_pw_pool)
do sic calculations on explicit orbitals
subroutine, public compute_matrix_vxc(qs_env, v_rspace, matrix_vxc)
compute matrix_vxc, defined via the potential created by qs_vxc_create ignores things like tau functi...
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kts, mu, flexible_electron_count)
Get the components of a MO set data structure.
superstucture that hold various representations of the density and keeps track of which ones are vali...
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
types for task lists
Exchange and Correlation functional calculations.
Definition xc.F:17
real(kind=dp) function, public xc_exc_calc(rho_r, rho_g, tau, xc_section, pw_pool)
calculates just the exchange and correlation energy (no vxc)
Definition xc.F:1456
subroutine, public xc_vxc_pw_create1(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, xc_section, pw_pool, compute_virial, virial_xc, exc_r)
Exchange and Correlation functional calculations. depending on the selected functional_routine calls ...
Definition xc.F:150
stores some data used in wavefunction fitting
Definition admm_types.F:120
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
keeps the information about the structure of a full matrix
represent a full matrix
type of a logger, at the moment it contains just a print level starting at which level it should be l...
stores some data used in construction of Kohn-Sham matrix
Definition hfx_types.F:510
Contains information about kpoints.
stores all the informations relevant to an mpi environment
contained for different pw related things
environment for the poisson solver
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Container for information about total charges on the grids.
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
keeps the density in various representations, keeping track of which ones are valid.