(git:374b731)
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-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief routines that build the Kohn-Sham matrix (i.e calculate the coulomb
10!> and xc parts
11!> \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
38 USE cp_fm_types, ONLY: cp_fm_create,&
47 USE cp_output_handling, ONLY: cp_p_file,&
51 USE dbcsr_api, ONLY: &
52 dbcsr_add, dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_dot, dbcsr_get_info, dbcsr_init_p, &
53 dbcsr_multiply, dbcsr_p_type, dbcsr_release_p, dbcsr_scale, dbcsr_scale_by_vector, &
54 dbcsr_set, dbcsr_type
57 USE hfx_types, ONLY: hfx_type
58 USE input_constants, ONLY: &
67 USE kinds, ONLY: default_string_length,&
68 dp
69 USE kpoint_types, ONLY: get_kpoint_info,&
80 USE ps_implicit_types, ONLY: mixed_bc,&
84 USE pw_env_types, ONLY: pw_env_get,&
86 USE pw_methods, ONLY: pw_axpy,&
87 pw_copy,&
90 pw_scale,&
97 USE pw_types, ONLY: pw_c1d_gs_type,&
106 USE qs_integrate_potential, ONLY: integrate_v_rspace,&
107 integrate_v_rspace_diagonal,&
108 integrate_v_rspace_one_center
109 USE qs_kind_types, ONLY: get_qs_kind_set,&
112 USE qs_ks_types, ONLY: get_ks_env,&
114 USE qs_mo_types, ONLY: get_mo_set,&
116 USE qs_rho_types, ONLY: qs_rho_get,&
119 USE virial_types, ONLY: virial_type
120 USE xc, ONLY: xc_exc_calc,&
122#include "./base/base_uses.f90"
123
124 IMPLICIT NONE
125
126 PRIVATE
127
128 LOGICAL, PARAMETER :: debug_this_module = .true.
129 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_utils'
130
134
135CONTAINS
136
137! **************************************************************************************************
138!> \brief do ROKS calculations yielding low spin states
139!> \param energy ...
140!> \param qs_env ...
141!> \param dft_control ...
142!> \param do_hfx ...
143!> \param just_energy ...
144!> \param calculate_forces ...
145!> \param auxbas_pw_pool ...
146! **************************************************************************************************
147 SUBROUTINE low_spin_roks(energy, qs_env, dft_control, do_hfx, just_energy, &
148 calculate_forces, auxbas_pw_pool)
149
150 TYPE(qs_energy_type), POINTER :: energy
151 TYPE(qs_environment_type), POINTER :: qs_env
152 TYPE(dft_control_type), POINTER :: dft_control
153 LOGICAL, INTENT(IN) :: do_hfx, just_energy, calculate_forces
154 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
155
156 CHARACTER(*), PARAMETER :: routinen = 'low_spin_roks'
157
158 INTEGER :: handle, irep, ispin, iterm, k, k_alpha, &
159 k_beta, n_rep, nelectron, nspin, nterms
160 INTEGER, DIMENSION(:), POINTER :: ivec
161 INTEGER, DIMENSION(:, :, :), POINTER :: occupations
162 LOGICAL :: compute_virial, in_range, &
163 uniform_occupation
164 REAL(kind=dp) :: ehfx, exc
165 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
166 REAL(kind=dp), DIMENSION(:), POINTER :: energy_scaling, rvec, scaling
167 TYPE(cell_type), POINTER :: cell
168 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_h, matrix_hfx, matrix_p, mdummy, &
169 mo_derivs, rho_ao
170 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p2
171 TYPE(dbcsr_type), POINTER :: dbcsr_deriv, fm_deriv, fm_scaled, &
172 mo_coeff
173 TYPE(hfx_type), DIMENSION(:, :), POINTER :: x_data
174 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
175 TYPE(mp_para_env_type), POINTER :: para_env
176 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
177 TYPE(pw_env_type), POINTER :: pw_env
178 TYPE(pw_pool_type), POINTER :: xc_pw_pool
179 TYPE(pw_r3d_rs_type) :: work_v_rspace
180 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc, vxc_tau
181 TYPE(qs_ks_env_type), POINTER :: ks_env
182 TYPE(qs_rho_type), POINTER :: rho
183 TYPE(section_vals_type), POINTER :: hfx_section, input, &
184 low_spin_roks_section, xc_section
185 TYPE(virial_type), POINTER :: virial
186
187 IF (.NOT. dft_control%low_spin_roks) RETURN
188
189 CALL timeset(routinen, handle)
190
191 NULLIFY (ks_env, rho_ao)
192
193 ! Test for not compatible options
194 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
195 CALL cp_abort(__location__, "GAPW/GAPW_XC are not compatible with low spin ROKS method.")
196 END IF
197 IF (dft_control%do_admm) THEN
198 CALL cp_abort(__location__, "ADMM not compatible with low spin ROKS method.")
199 END IF
200 IF (dft_control%do_admm) THEN
201 IF (qs_env%admm_env%aux_exch_func /= do_admm_aux_exch_func_none) THEN
202 CALL cp_abort(__location__, "ADMM with XC correction functional "// &
203 "not compatible with low spin ROKS method.")
204 END IF
205 END IF
206 IF (dft_control%qs_control%semi_empirical .OR. dft_control%qs_control%dftb .OR. &
207 dft_control%qs_control%xtb) THEN
208 CALL cp_abort(__location__, "SE/xTB/DFTB are not compatible with low spin ROKS method.")
209 END IF
210
211 CALL get_qs_env(qs_env, &
212 ks_env=ks_env, &
213 mo_derivs=mo_derivs, &
214 mos=mo_array, &
215 rho=rho, &
216 pw_env=pw_env, &
217 input=input, &
218 cell=cell, &
219 virial=virial)
220
221 CALL qs_rho_get(rho, rho_ao=rho_ao)
222
223 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
224 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
225 hfx_section => section_vals_get_subs_vals(input, "DFT%XC%HF")
226
227 ! some assumptions need to be checked
228 ! we have two spins
229 cpassert(SIZE(mo_array, 1) == 2)
230 nspin = 2
231 ! we want uniform occupations
232 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
233 cpassert(uniform_occupation)
234 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff, uniform_occupation=uniform_occupation)
235 cpassert(uniform_occupation)
236 IF (do_hfx .AND. calculate_forces .AND. compute_virial) THEN
237 CALL cp_abort(__location__, "ROKS virial with HFX not available.")
238 END IF
239
240 NULLIFY (dbcsr_deriv)
241 CALL dbcsr_init_p(dbcsr_deriv)
242 CALL dbcsr_copy(dbcsr_deriv, mo_derivs(1)%matrix)
243 CALL dbcsr_set(dbcsr_deriv, 0.0_dp)
244
245 ! basic info
246 CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
247 CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_alpha)
248 CALL get_mo_set(mo_set=mo_array(2), mo_coeff_b=mo_coeff)
249 CALL dbcsr_get_info(mo_coeff, nfullcols_total=k_beta)
250
251 ! read the input
252 low_spin_roks_section => section_vals_get_subs_vals(input, "DFT%LOW_SPIN_ROKS")
253
254 CALL section_vals_val_get(low_spin_roks_section, "ENERGY_SCALING", r_vals=rvec)
255 nterms = SIZE(rvec)
256 ALLOCATE (energy_scaling(nterms))
257 energy_scaling = rvec !? just wondering, should this add up to 1, in which case we should cpp?
258
259 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", n_rep_val=n_rep)
260 cpassert(n_rep == nterms)
261 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=1, i_vals=ivec)
262 nelectron = SIZE(ivec)
263 cpassert(nelectron == k_alpha - k_beta)
264 ALLOCATE (occupations(2, nelectron, nterms))
265 occupations = 0
266 DO iterm = 1, nterms
267 CALL section_vals_val_get(low_spin_roks_section, "SPIN_CONFIGURATION", i_rep_val=iterm, i_vals=ivec)
268 cpassert(nelectron == SIZE(ivec))
269 in_range = all(ivec >= 1) .AND. all(ivec <= 2)
270 cpassert(in_range)
271 DO k = 1, nelectron
272 occupations(ivec(k), k, iterm) = 1
273 END DO
274 END DO
275
276 ! set up general data structures
277 ! density matrices, kohn-sham matrices
278
279 NULLIFY (matrix_p)
280 CALL dbcsr_allocate_matrix_set(matrix_p, nspin)
281 DO ispin = 1, nspin
282 ALLOCATE (matrix_p(ispin)%matrix)
283 CALL dbcsr_copy(matrix_p(ispin)%matrix, rho_ao(1)%matrix, &
284 name="density matrix low spin roks")
285 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
286 END DO
287
288 NULLIFY (matrix_h)
289 CALL dbcsr_allocate_matrix_set(matrix_h, nspin)
290 DO ispin = 1, nspin
291 ALLOCATE (matrix_h(ispin)%matrix)
292 CALL dbcsr_copy(matrix_h(ispin)%matrix, rho_ao(1)%matrix, &
293 name="KS matrix low spin roks")
294 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
295 END DO
296
297 IF (do_hfx) THEN
298 NULLIFY (matrix_hfx)
299 CALL dbcsr_allocate_matrix_set(matrix_hfx, nspin)
300 DO ispin = 1, nspin
301 ALLOCATE (matrix_hfx(ispin)%matrix)
302 CALL dbcsr_copy(matrix_hfx(ispin)%matrix, rho_ao(1)%matrix, &
303 name="HFX matrix low spin roks")
304 END DO
305 END IF
306
307 ! grids in real and g space for rho and vxc
308 ! tau functionals are not supported
309 NULLIFY (tau, vxc_tau, vxc)
310 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
311
312 ALLOCATE (rho_r(nspin))
313 ALLOCATE (rho_g(nspin))
314 DO ispin = 1, nspin
315 CALL auxbas_pw_pool%create_pw(rho_r(ispin))
316 CALL auxbas_pw_pool%create_pw(rho_g(ispin))
317 END DO
318 CALL auxbas_pw_pool%create_pw(work_v_rspace)
319
320 ! get mo matrices needed to construct the density matrices
321 ! we will base all on the alpha spin matrix, obviously possible in ROKS
322 CALL get_mo_set(mo_set=mo_array(1), mo_coeff_b=mo_coeff)
323 NULLIFY (fm_scaled, fm_deriv)
324 CALL dbcsr_init_p(fm_scaled)
325 CALL dbcsr_init_p(fm_deriv)
326 CALL dbcsr_copy(fm_scaled, mo_coeff)
327 CALL dbcsr_copy(fm_deriv, mo_coeff)
328
329 ALLOCATE (scaling(k_alpha))
330
331 ! for each term, add it with the given scaling factor to the energy, and compute the required derivatives
332 DO iterm = 1, nterms
333
334 DO ispin = 1, nspin
335 ! compute the proper density matrices with the required occupations
336 CALL dbcsr_set(matrix_p(ispin)%matrix, 0.0_dp)
337 scaling = 1.0_dp
338 scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
339 CALL dbcsr_copy(fm_scaled, mo_coeff)
340 CALL dbcsr_scale_by_vector(fm_scaled, scaling, side='right')
341 CALL dbcsr_multiply('n', 't', 1.0_dp, mo_coeff, fm_scaled, &
342 0.0_dp, matrix_p(ispin)%matrix, retain_sparsity=.true.)
343 ! compute the densities on the grid
344 CALL calculate_rho_elec(matrix_p=matrix_p(ispin)%matrix, &
345 rho=rho_r(ispin), rho_gspace=rho_g(ispin), &
346 ks_env=ks_env)
347 END DO
348
349 ! compute the exchange energies / potential if needed
350 IF (just_energy) THEN
351 exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
352 pw_pool=xc_pw_pool)
353 ELSE
354 cpassert(.NOT. compute_virial)
355 CALL xc_vxc_pw_create1(vxc_rho=vxc, rho_r=rho_r, &
356 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
357 pw_pool=xc_pw_pool, compute_virial=.false., virial_xc=virial_xc_tmp)
358 END IF
359
360 energy%exc = energy%exc + energy_scaling(iterm)*exc
361
362 IF (do_hfx) THEN
363 ! Add Hartree-Fock contribution
364 DO ispin = 1, nspin
365 CALL dbcsr_set(matrix_hfx(ispin)%matrix, 0.0_dp)
366 END DO
367 ehfx = energy%ex
368 CALL tddft_hfx_matrix(matrix_hfx, matrix_p, qs_env, &
369 recalc_integrals=.false., update_energy=.true.)
370 energy%ex = ehfx + energy_scaling(iterm)*energy%ex
371 END IF
372
373 ! add the corresponding derivatives to the MO derivatives
374 IF (.NOT. just_energy) THEN
375 ! get the potential in matrix form
376 DO ispin = 1, nspin
377 CALL dbcsr_set(matrix_h(ispin)%matrix, 0.0_dp)
378 ! use a work_v_rspace
379 CALL pw_axpy(vxc(ispin), work_v_rspace, energy_scaling(iterm)*vxc(ispin)%pw_grid%dvol, 0.0_dp)
380 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=matrix_p(ispin), hmat=matrix_h(ispin), &
381 qs_env=qs_env, calculate_forces=calculate_forces)
382 CALL auxbas_pw_pool%give_back_pw(vxc(ispin))
383 END DO
384 DEALLOCATE (vxc)
385
386 IF (do_hfx) THEN
387 ! add HFX contribution
388 DO ispin = 1, nspin
389 CALL dbcsr_add(matrix_h(ispin)%matrix, matrix_hfx(ispin)%matrix, &
390 1.0_dp, energy_scaling(iterm))
391 END DO
392 IF (calculate_forces) THEN
393 CALL get_qs_env(qs_env, x_data=x_data, para_env=para_env)
394 IF (x_data(1, 1)%n_rep_hf /= 1) THEN
395 CALL cp_abort(__location__, "Multiple HFX section forces not compatible "// &
396 "with low spin ROKS method.")
397 END IF
398 IF (x_data(1, 1)%do_hfx_ri) THEN
399 CALL cp_abort(__location__, "HFX_RI forces not compatible with low spin ROKS method.")
400 ELSE
401 irep = 1
402 NULLIFY (mdummy)
403 matrix_p2(1:nspin, 1:1) => matrix_p(1:nspin)
404 CALL derivatives_four_center(qs_env, matrix_p2, mdummy, hfx_section, para_env, &
405 irep, compute_virial, &
406 adiabatic_rescale_factor=energy_scaling(iterm))
407 END IF
408 END IF
409
410 END IF
411
412 ! add this to the mo_derivs, again based on the alpha mo_coeff
413 DO ispin = 1, nspin
414 CALL dbcsr_multiply('n', 'n', 1.0_dp, matrix_h(ispin)%matrix, mo_coeff, &
415 0.0_dp, dbcsr_deriv, last_column=k_alpha)
416
417 scaling = 1.0_dp
418 scaling(k_alpha - nelectron + 1:k_alpha) = occupations(ispin, :, iterm)
419 CALL dbcsr_scale_by_vector(dbcsr_deriv, scaling, side='right')
420 CALL dbcsr_add(mo_derivs(1)%matrix, dbcsr_deriv, 1.0_dp, 1.0_dp)
421 END DO
422
423 END IF
424
425 END DO
426
427 ! release allocated memory
428 DO ispin = 1, nspin
429 CALL auxbas_pw_pool%give_back_pw(rho_r(ispin))
430 CALL auxbas_pw_pool%give_back_pw(rho_g(ispin))
431 END DO
432 DEALLOCATE (rho_r, rho_g)
433 CALL dbcsr_deallocate_matrix_set(matrix_p)
434 CALL dbcsr_deallocate_matrix_set(matrix_h)
435 IF (do_hfx) THEN
436 CALL dbcsr_deallocate_matrix_set(matrix_hfx)
437 END IF
438
439 CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
440
441 CALL dbcsr_release_p(fm_deriv)
442 CALL dbcsr_release_p(fm_scaled)
443
444 DEALLOCATE (occupations)
445 DEALLOCATE (energy_scaling)
446 DEALLOCATE (scaling)
447
448 CALL dbcsr_release_p(dbcsr_deriv)
449
450 CALL timestop(handle)
451
452 END SUBROUTINE low_spin_roks
453! **************************************************************************************************
454!> \brief do sic calculations on explicit orbitals
455!> \param energy ...
456!> \param qs_env ...
457!> \param dft_control ...
458!> \param poisson_env ...
459!> \param just_energy ...
460!> \param calculate_forces ...
461!> \param auxbas_pw_pool ...
462! **************************************************************************************************
463 SUBROUTINE sic_explicit_orbitals(energy, qs_env, dft_control, poisson_env, just_energy, &
464 calculate_forces, auxbas_pw_pool)
465
466 TYPE(qs_energy_type), POINTER :: energy
467 TYPE(qs_environment_type), POINTER :: qs_env
468 TYPE(dft_control_type), POINTER :: dft_control
469 TYPE(pw_poisson_type), POINTER :: poisson_env
470 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
471 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
472
473 CHARACTER(*), PARAMETER :: routinen = 'sic_explicit_orbitals'
474
475 INTEGER :: handle, i, iorb, k_alpha, k_beta, norb
476 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: sic_orbital_list
477 LOGICAL :: compute_virial, uniform_occupation
478 REAL(kind=dp) :: ener, exc
479 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
480 TYPE(cell_type), POINTER :: cell
481 TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
482 TYPE(cp_fm_type) :: matrix_hv, matrix_v
483 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: mo_derivs_local
484 TYPE(cp_fm_type), POINTER :: mo_coeff
485 TYPE(dbcsr_p_type) :: orb_density_matrix_p, orb_h_p
486 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mo_derivs, rho_ao, tmp_dbcsr
487 TYPE(dbcsr_type), POINTER :: orb_density_matrix, orb_h
488 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
489 TYPE(pw_c1d_gs_type) :: work_v_gspace
490 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
491 TYPE(pw_c1d_gs_type), TARGET :: orb_rho_g, tmp_g
492 TYPE(pw_env_type), POINTER :: pw_env
493 TYPE(pw_pool_type), POINTER :: xc_pw_pool
494 TYPE(pw_r3d_rs_type) :: work_v_rspace
495 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau, vxc, vxc_tau
496 TYPE(pw_r3d_rs_type), TARGET :: orb_rho_r, tmp_r
497 TYPE(qs_ks_env_type), POINTER :: ks_env
498 TYPE(qs_rho_type), POINTER :: rho
499 TYPE(section_vals_type), POINTER :: input, xc_section
500 TYPE(virial_type), POINTER :: virial
501
502 IF (dft_control%sic_method_id .NE. sic_eo) RETURN
503
504 CALL timeset(routinen, handle)
505
506 NULLIFY (tau, vxc_tau, mo_derivs, ks_env, rho_ao)
507
508 ! generate the lists of orbitals that need sic treatment
509 CALL get_qs_env(qs_env, &
510 ks_env=ks_env, &
511 mo_derivs=mo_derivs, &
512 mos=mo_array, &
513 rho=rho, &
514 pw_env=pw_env, &
515 input=input, &
516 cell=cell, &
517 virial=virial)
518
519 CALL qs_rho_get(rho, rho_ao=rho_ao)
520
521 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
522 xc_section => section_vals_get_subs_vals(input, "DFT%XC")
523
524 DO i = 1, SIZE(mo_array) !fm->dbcsr
525 IF (mo_array(i)%use_mo_coeff_b) THEN !fm->dbcsr
526 CALL copy_dbcsr_to_fm(mo_array(i)%mo_coeff_b, &
527 mo_array(i)%mo_coeff) !fm->dbcsr
528 END IF !fm->dbcsr
529 END DO !fm->dbcsr
530
531 CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool)
532
533 ! we have two spins
534 cpassert(SIZE(mo_array, 1) == 2)
535 ! we want uniform occupations
536 CALL get_mo_set(mo_set=mo_array(1), uniform_occupation=uniform_occupation)
537 cpassert(uniform_occupation)
538 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff, uniform_occupation=uniform_occupation)
539 cpassert(uniform_occupation)
540
541 NULLIFY (tmp_dbcsr)
542 CALL dbcsr_allocate_matrix_set(tmp_dbcsr, SIZE(mo_derivs, 1))
543 DO i = 1, SIZE(mo_derivs, 1) !fm->dbcsr
544 !
545 NULLIFY (tmp_dbcsr(i)%matrix)
546 CALL dbcsr_init_p(tmp_dbcsr(i)%matrix)
547 CALL dbcsr_copy(tmp_dbcsr(i)%matrix, mo_derivs(i)%matrix)
548 CALL dbcsr_set(tmp_dbcsr(i)%matrix, 0.0_dp)
549 END DO !fm->dbcsr
550
551 k_alpha = 0; k_beta = 0
552 SELECT CASE (dft_control%sic_list_id)
553 CASE (sic_list_all)
554
555 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
556 CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
557
558 IF (SIZE(mo_array, 1) > 1) THEN
559 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
560 CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
561 END IF
562
563 norb = k_alpha + k_beta
564 ALLOCATE (sic_orbital_list(3, norb))
565
566 iorb = 0
567 DO i = 1, k_alpha
568 iorb = iorb + 1
569 sic_orbital_list(1, iorb) = 1
570 sic_orbital_list(2, iorb) = i
571 sic_orbital_list(3, iorb) = 1
572 END DO
573 DO i = 1, k_beta
574 iorb = iorb + 1
575 sic_orbital_list(1, iorb) = 2
576 sic_orbital_list(2, iorb) = i
577 IF (SIZE(mo_derivs, 1) == 1) THEN
578 sic_orbital_list(3, iorb) = 1
579 ELSE
580 sic_orbital_list(3, iorb) = 2
581 END IF
582 END DO
583
584 CASE (sic_list_unpaired)
585 ! we have two spins
586 cpassert(SIZE(mo_array, 1) == 2)
587 ! we have them restricted
588 cpassert(SIZE(mo_derivs, 1) == 1)
589 cpassert(dft_control%restricted)
590
591 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
592 CALL cp_fm_get_info(mo_coeff, ncol_global=k_alpha)
593
594 CALL get_mo_set(mo_set=mo_array(2), mo_coeff=mo_coeff)
595 CALL cp_fm_get_info(mo_coeff, ncol_global=k_beta)
596
597 norb = k_alpha - k_beta
598 ALLOCATE (sic_orbital_list(3, norb))
599
600 iorb = 0
601 DO i = k_beta + 1, k_alpha
602 iorb = iorb + 1
603 sic_orbital_list(1, iorb) = 1
604 sic_orbital_list(2, iorb) = i
605 ! we are guaranteed to be restricted
606 sic_orbital_list(3, iorb) = 1
607 END DO
608
609 CASE DEFAULT
610 cpabort("")
611 END SELECT
612
613 ! data needed for each of the orbs
614 CALL auxbas_pw_pool%create_pw(orb_rho_r)
615 CALL auxbas_pw_pool%create_pw(tmp_r)
616 CALL auxbas_pw_pool%create_pw(orb_rho_g)
617 CALL auxbas_pw_pool%create_pw(tmp_g)
618 CALL auxbas_pw_pool%create_pw(work_v_gspace)
619 CALL auxbas_pw_pool%create_pw(work_v_rspace)
620
621 ALLOCATE (orb_density_matrix)
622 CALL dbcsr_copy(orb_density_matrix, rho_ao(1)%matrix, &
623 name="orb_density_matrix")
624 CALL dbcsr_set(orb_density_matrix, 0.0_dp)
625 orb_density_matrix_p%matrix => orb_density_matrix
626
627 ALLOCATE (orb_h)
628 CALL dbcsr_copy(orb_h, rho_ao(1)%matrix, &
629 name="orb_density_matrix")
630 CALL dbcsr_set(orb_h, 0.0_dp)
631 orb_h_p%matrix => orb_h
632
633 CALL get_mo_set(mo_set=mo_array(1), mo_coeff=mo_coeff)
634
635 CALL cp_fm_struct_create(fm_struct_tmp, ncol_global=1, &
636 template_fmstruct=mo_coeff%matrix_struct)
637 CALL cp_fm_create(matrix_v, fm_struct_tmp, name="matrix_v")
638 CALL cp_fm_create(matrix_hv, fm_struct_tmp, name="matrix_hv")
639 CALL cp_fm_struct_release(fm_struct_tmp)
640
641 ALLOCATE (mo_derivs_local(SIZE(mo_array, 1)))
642 DO i = 1, SIZE(mo_array, 1)
643 CALL get_mo_set(mo_set=mo_array(i), mo_coeff=mo_coeff)
644 CALL cp_fm_create(mo_derivs_local(i), mo_coeff%matrix_struct)
645 END DO
646
647 ALLOCATE (rho_r(2))
648 rho_r(1) = orb_rho_r
649 rho_r(2) = tmp_r
650 CALL pw_zero(tmp_r)
651
652 ALLOCATE (rho_g(2))
653 rho_g(1) = orb_rho_g
654 rho_g(2) = tmp_g
655 CALL pw_zero(tmp_g)
656
657 NULLIFY (vxc)
658 ! now apply to SIC correction to each selected orbital
659 DO iorb = 1, norb
660 ! extract the proper orbital from the mo_coeff
661 CALL get_mo_set(mo_set=mo_array(sic_orbital_list(1, iorb)), mo_coeff=mo_coeff)
662 CALL cp_fm_to_fm(mo_coeff, matrix_v, 1, sic_orbital_list(2, iorb), 1)
663
664 ! construct the density matrix and the corresponding density
665 CALL dbcsr_set(orb_density_matrix, 0.0_dp)
666 CALL cp_dbcsr_plus_fm_fm_t(orb_density_matrix, matrix_v=matrix_v, ncol=1, &
667 alpha=1.0_dp)
668
669 CALL calculate_rho_elec(matrix_p=orb_density_matrix, &
670 rho=orb_rho_r, rho_gspace=orb_rho_g, &
671 ks_env=ks_env)
672
673 ! compute the energy functional for this orbital and its derivative
674
675 CALL pw_poisson_solve(poisson_env, orb_rho_g, ener, work_v_gspace)
676 ! no PBC correction is done here, see "calc_v_sic_rspace" for SIC methods
677 ! with PBC aware corrections
678 energy%hartree = energy%hartree - dft_control%sic_scaling_a*ener
679 IF (.NOT. just_energy) THEN
680 CALL pw_transfer(work_v_gspace, work_v_rspace)
681 CALL pw_scale(work_v_rspace, -dft_control%sic_scaling_a*work_v_rspace%pw_grid%dvol)
682 CALL dbcsr_set(orb_h, 0.0_dp)
683 END IF
684
685 IF (just_energy) THEN
686 exc = xc_exc_calc(rho_r=rho_r, rho_g=rho_g, tau=tau, xc_section=xc_section, &
687 pw_pool=xc_pw_pool)
688 ELSE
689 cpassert(.NOT. compute_virial)
690 CALL xc_vxc_pw_create1(vxc_rho=vxc, rho_r=rho_r, &
691 rho_g=rho_g, tau=tau, vxc_tau=vxc_tau, exc=exc, xc_section=xc_section, &
692 pw_pool=xc_pw_pool, compute_virial=compute_virial, virial_xc=virial_xc_tmp)
693 ! add to the existing work_v_rspace
694 CALL pw_axpy(vxc(1), work_v_rspace, -dft_control%sic_scaling_b*vxc(1)%pw_grid%dvol)
695 END IF
696 energy%exc = energy%exc - dft_control%sic_scaling_b*exc
697
698 IF (.NOT. just_energy) THEN
699 ! note, orb_h (which is being pointed to with orb_h_p) is zeroed above
700 CALL integrate_v_rspace(v_rspace=work_v_rspace, pmat=orb_density_matrix_p, hmat=orb_h_p, &
701 qs_env=qs_env, calculate_forces=calculate_forces)
702
703 ! add this to the mo_derivs
704 CALL cp_dbcsr_sm_fm_multiply(orb_h, matrix_v, matrix_hv, 1)
705 ! silly trick, copy to an array of the right size and add to mo_derivs
706 CALL cp_fm_set_all(mo_derivs_local(sic_orbital_list(3, iorb)), 0.0_dp)
707 CALL cp_fm_to_fm(matrix_hv, mo_derivs_local(sic_orbital_list(3, iorb)), 1, 1, sic_orbital_list(2, iorb))
708 CALL copy_fm_to_dbcsr(mo_derivs_local(sic_orbital_list(3, iorb)), &
709 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix)
710 CALL dbcsr_add(mo_derivs(sic_orbital_list(3, iorb))%matrix, &
711 tmp_dbcsr(sic_orbital_list(3, iorb))%matrix, 1.0_dp, 1.0_dp)
712 !
713 ! need to deallocate vxc
714 CALL xc_pw_pool%give_back_pw(vxc(1))
715 CALL xc_pw_pool%give_back_pw(vxc(2))
716 DEALLOCATE (vxc)
717
718 END IF
719
720 END DO
721
722 CALL auxbas_pw_pool%give_back_pw(orb_rho_r)
723 CALL auxbas_pw_pool%give_back_pw(tmp_r)
724 CALL auxbas_pw_pool%give_back_pw(orb_rho_g)
725 CALL auxbas_pw_pool%give_back_pw(tmp_g)
726 CALL auxbas_pw_pool%give_back_pw(work_v_gspace)
727 CALL auxbas_pw_pool%give_back_pw(work_v_rspace)
728
729 CALL dbcsr_deallocate_matrix(orb_density_matrix)
730 CALL dbcsr_deallocate_matrix(orb_h)
731 CALL cp_fm_release(matrix_v)
732 CALL cp_fm_release(matrix_hv)
733 CALL cp_fm_release(mo_derivs_local)
734 DEALLOCATE (rho_r)
735 DEALLOCATE (rho_g)
736
737 CALL dbcsr_deallocate_matrix_set(tmp_dbcsr) !fm->dbcsr
738
739 CALL timestop(handle)
740
741 END SUBROUTINE sic_explicit_orbitals
742
743! **************************************************************************************************
744!> \brief do sic calculations on the spin density
745!> \param v_sic_rspace ...
746!> \param energy ...
747!> \param qs_env ...
748!> \param dft_control ...
749!> \param rho ...
750!> \param poisson_env ...
751!> \param just_energy ...
752!> \param calculate_forces ...
753!> \param auxbas_pw_pool ...
754! **************************************************************************************************
755 SUBROUTINE calc_v_sic_rspace(v_sic_rspace, energy, &
756 qs_env, dft_control, rho, poisson_env, just_energy, &
757 calculate_forces, auxbas_pw_pool)
758
759 TYPE(pw_r3d_rs_type), POINTER :: v_sic_rspace
760 TYPE(qs_energy_type), POINTER :: energy
761 TYPE(qs_environment_type), POINTER :: qs_env
762 TYPE(dft_control_type), POINTER :: dft_control
763 TYPE(qs_rho_type), POINTER :: rho
764 TYPE(pw_poisson_type), POINTER :: poisson_env
765 LOGICAL, INTENT(IN) :: just_energy, calculate_forces
766 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
767
768 INTEGER :: i, nelec, nelec_a, nelec_b, nforce
769 REAL(kind=dp) :: ener, full_scaling, scaling
770 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: store_forces
771 TYPE(mo_set_type), DIMENSION(:), POINTER :: mo_array
772 TYPE(pw_c1d_gs_type) :: work_rho, work_v
773 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
774 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
775
776 NULLIFY (mo_array, rho_g)
777
778 IF (dft_control%sic_method_id == sic_none) RETURN
779 IF (dft_control%sic_method_id == sic_eo) RETURN
780
781 IF (dft_control%qs_control%gapw) &
782 cpabort("sic and GAPW not yet compatible")
783
784 ! OK, right now we like two spins to do sic, could be relaxed for AD
785 cpassert(dft_control%nspins == 2)
786
787 CALL auxbas_pw_pool%create_pw(work_rho)
788 CALL auxbas_pw_pool%create_pw(work_v)
789
790 CALL qs_rho_get(rho, rho_g=rho_g)
791
792 ! Hartree sic corrections
793 SELECT CASE (dft_control%sic_method_id)
795 CALL pw_copy(rho_g(1), work_rho)
796 CALL pw_axpy(rho_g(2), work_rho, alpha=-1._dp)
797 CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
798 CASE (sic_ad)
799 ! find out how many elecs we have
800 CALL get_qs_env(qs_env, mos=mo_array)
801 CALL get_mo_set(mo_set=mo_array(1), nelectron=nelec_a)
802 CALL get_mo_set(mo_set=mo_array(2), nelectron=nelec_b)
803 nelec = nelec_a + nelec_b
804 CALL pw_copy(rho_g(1), work_rho)
805 CALL pw_axpy(rho_g(2), work_rho)
806 scaling = 1.0_dp/real(nelec, kind=dp)
807 CALL pw_scale(work_rho, scaling)
808 CALL pw_poisson_solve(poisson_env, work_rho, ener, work_v)
809 CASE DEFAULT
810 cpabort("Unknown sic method id")
811 END SELECT
812
813 ! Correct for DDAP charges (if any)
814 ! storing whatever force might be there from previous decoupling
815 IF (calculate_forces) THEN
816 CALL get_qs_env(qs_env=qs_env, force=force)
817 nforce = 0
818 DO i = 1, SIZE(force)
819 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
820 END DO
821 ALLOCATE (store_forces(3, nforce))
822 nforce = 0
823 DO i = 1, SIZE(force)
824 store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2)) = force(i)%ch_pulay(:, :)
825 force(i)%ch_pulay(:, :) = 0.0_dp
826 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
827 END DO
828 END IF
829
830 CALL cp_ddapc_apply_cd(qs_env, &
831 work_rho, &
832 ener, &
833 v_hartree_gspace=work_v, &
834 calculate_forces=calculate_forces, &
835 itype_of_density="SPIN")
836
837 SELECT CASE (dft_control%sic_method_id)
839 full_scaling = -dft_control%sic_scaling_a
840 CASE (sic_ad)
841 full_scaling = -dft_control%sic_scaling_a*nelec
842 CASE DEFAULT
843 cpabort("Unknown sic method id")
844 END SELECT
845 energy%hartree = energy%hartree + full_scaling*ener
846
847 ! add scaled forces, restoring the old
848 IF (calculate_forces) THEN
849 nforce = 0
850 DO i = 1, SIZE(force)
851 force(i)%ch_pulay(:, :) = force(i)%ch_pulay(:, :)*full_scaling + &
852 store_forces(1:3, nforce + 1:nforce + SIZE(force(i)%ch_pulay, 2))
853 nforce = nforce + SIZE(force(i)%ch_pulay, 2)
854 END DO
855 END IF
856
857 IF (.NOT. just_energy) THEN
858 ALLOCATE (v_sic_rspace)
859 CALL auxbas_pw_pool%create_pw(v_sic_rspace)
860 CALL pw_transfer(work_v, v_sic_rspace)
861 ! also take into account the scaling (in addition to the volume element)
862 CALL pw_scale(v_sic_rspace, &
863 dft_control%sic_scaling_a*v_sic_rspace%pw_grid%dvol)
864 END IF
865
866 CALL auxbas_pw_pool%give_back_pw(work_rho)
867 CALL auxbas_pw_pool%give_back_pw(work_v)
868
869 END SUBROUTINE calc_v_sic_rspace
870
871! **************************************************************************************************
872!> \brief ...
873!> \param qs_env ...
874!> \param rho ...
875! **************************************************************************************************
876 SUBROUTINE print_densities(qs_env, rho)
877 TYPE(qs_environment_type), POINTER :: qs_env
878 TYPE(qs_rho_type), POINTER :: rho
879
880 INTEGER :: img, ispin, n_electrons, output_unit
881 REAL(dp) :: tot1_h, tot1_s, tot_rho_r, trace, &
882 trace_tmp
883 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r_arr
884 TYPE(cell_type), POINTER :: cell
885 TYPE(cp_logger_type), POINTER :: logger
886 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s, rho_ao
887 TYPE(dft_control_type), POINTER :: dft_control
888 TYPE(qs_charges_type), POINTER :: qs_charges
889 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
890 TYPE(section_vals_type), POINTER :: input, scf_section
891
892 NULLIFY (qs_charges, qs_kind_set, cell, input, logger, scf_section, matrix_s, &
893 dft_control, tot_rho_r_arr, rho_ao)
894
895 logger => cp_get_default_logger()
896
897 CALL get_qs_env(qs_env, &
898 qs_kind_set=qs_kind_set, &
899 cell=cell, qs_charges=qs_charges, &
900 input=input, &
901 matrix_s_kp=matrix_s, &
902 dft_control=dft_control)
903
904 CALL get_qs_kind_set(qs_kind_set, nelectron=n_electrons)
905
906 scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
907 output_unit = cp_print_key_unit_nr(logger, scf_section, "PRINT%TOTAL_DENSITIES", &
908 extension=".scfLog")
909
910 CALL qs_rho_get(rho, tot_rho_r=tot_rho_r_arr, rho_ao_kp=rho_ao)
911 n_electrons = n_electrons - dft_control%charge
912 tot_rho_r = accurate_sum(tot_rho_r_arr)
913
914 trace = 0
915 IF (btest(cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES"), cp_p_file)) THEN
916 DO ispin = 1, dft_control%nspins
917 DO img = 1, dft_control%nimages
918 CALL dbcsr_dot(rho_ao(ispin, img)%matrix, matrix_s(1, img)%matrix, trace_tmp)
919 trace = trace + trace_tmp
920 END DO
921 END DO
922 END IF
923
924 IF (output_unit > 0) THEN
925 WRITE (unit=output_unit, fmt="(/,T3,A,T41,F20.10)") "Trace(PS):", trace
926 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
927 "Electronic density on regular grids: ", &
928 tot_rho_r, &
929 tot_rho_r + &
930 REAL(n_electrons, dp), &
931 "Core density on regular grids:", &
932 qs_charges%total_rho_core_rspace, &
933 qs_charges%total_rho_core_rspace - REAL(n_electrons + dft_control%charge, dp)
934 END IF
935 IF (dft_control%qs_control%gapw) THEN
936 tot1_h = qs_charges%total_rho1_hard(1)
937 tot1_s = qs_charges%total_rho1_soft(1)
938 DO ispin = 2, dft_control%nspins
939 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
940 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
941 END DO
942 IF (output_unit > 0) THEN
943 WRITE (unit=output_unit, fmt="((T3,A,T41,2F20.10))") &
944 "Hard and soft densities (Lebedev):", &
945 tot1_h, tot1_s
946 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
947 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
948 tot_rho_r + tot1_h - tot1_s, &
949 "Total charge density (r-space): ", &
950 tot_rho_r + tot1_h - tot1_s &
951 + qs_charges%total_rho_core_rspace, &
952 "Total Rho_soft + Rho0_soft (g-space):", &
953 qs_charges%total_rho_gspace
954 END IF
955 qs_charges%background = tot_rho_r + tot1_h - tot1_s + &
956 qs_charges%total_rho_core_rspace
957 ELSE IF (dft_control%qs_control%gapw_xc) THEN
958 tot1_h = qs_charges%total_rho1_hard(1)
959 tot1_s = qs_charges%total_rho1_soft(1)
960 DO ispin = 2, dft_control%nspins
961 tot1_h = tot1_h + qs_charges%total_rho1_hard(ispin)
962 tot1_s = tot1_s + qs_charges%total_rho1_soft(ispin)
963 END DO
964 IF (output_unit > 0) THEN
965 WRITE (unit=output_unit, fmt="(/,(T3,A,T41,2F20.10))") &
966 "Hard and soft densities (Lebedev):", &
967 tot1_h, tot1_s
968 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
969 "Total Rho_soft + Rho1_hard - Rho1_soft (r-space): ", &
970 accurate_sum(tot_rho_r_arr) + tot1_h - tot1_s
971 END IF
972 qs_charges%background = tot_rho_r + &
973 qs_charges%total_rho_core_rspace
974 ELSE
975 IF (output_unit > 0) THEN
976 WRITE (unit=output_unit, fmt="(T3,A,T41,F20.10)") &
977 "Total charge density on r-space grids: ", &
978 tot_rho_r + &
979 qs_charges%total_rho_core_rspace, &
980 "Total charge density g-space grids: ", &
981 qs_charges%total_rho_gspace
982 END IF
983 qs_charges%background = tot_rho_r + &
984 qs_charges%total_rho_core_rspace
985 END IF
986 IF (output_unit > 0) WRITE (unit=output_unit, fmt="()")
987 qs_charges%background = qs_charges%background/cell%deth
988
989 CALL cp_print_key_finished_output(output_unit, logger, scf_section, &
990 "PRINT%TOTAL_DENSITIES")
991
992 END SUBROUTINE print_densities
993
994! **************************************************************************************************
995!> \brief Print detailed energies
996!>
997!> \param qs_env ...
998!> \param dft_control ...
999!> \param input ...
1000!> \param energy ...
1001!> \param mulliken_order_p ...
1002!> \par History
1003!> refactoring 04.03.2011 [MI]
1004!> \author
1005! **************************************************************************************************
1006 SUBROUTINE print_detailed_energy(qs_env, dft_control, input, energy, mulliken_order_p)
1007
1008 TYPE(qs_environment_type), POINTER :: qs_env
1009 TYPE(dft_control_type), POINTER :: dft_control
1010 TYPE(section_vals_type), POINTER :: input
1011 TYPE(qs_energy_type), POINTER :: energy
1012 REAL(kind=dp), INTENT(IN) :: mulliken_order_p
1013
1014 INTEGER :: bc, n, output_unit, psolver
1015 REAL(kind=dp) :: ddapc_order_p, implicit_ps_ehartree, &
1016 s2_order_p
1017 TYPE(cp_logger_type), POINTER :: logger
1018 TYPE(pw_env_type), POINTER :: pw_env
1019
1020 logger => cp_get_default_logger()
1021
1022 NULLIFY (pw_env)
1023 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1024 psolver = pw_env%poisson_env%parameters%solver
1025
1026 output_unit = cp_print_key_unit_nr(logger, input, "DFT%SCF%PRINT%DETAILED_ENERGY", &
1027 extension=".scfLog")
1028 IF (output_unit > 0) THEN
1029 IF (dft_control%do_admm) THEN
1030 WRITE (unit=output_unit, fmt="((T3,A,T60,F20.10))") &
1031 "Wfn fit exchange-correlation energy: ", energy%exc_aux_fit
1032 IF (dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1033 WRITE (unit=output_unit, fmt="((T3,A,T60,F20.10))") &
1034 "Wfn fit soft/hard atomic rho1 Exc contribution: ", energy%exc1_aux_fit
1035 END IF
1036 END IF
1037 IF (dft_control%do_admm) THEN
1038 IF (psolver .EQ. pw_poisson_implicit) THEN
1039 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1040 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1041 SELECT CASE (bc)
1043 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1044 "Core Hamiltonian energy: ", energy%core, &
1045 "Hartree energy: ", implicit_ps_ehartree, &
1046 "Electric enthalpy: ", energy%hartree, &
1047 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1048 CASE (periodic_bc, neumann_bc)
1049 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1050 "Core Hamiltonian energy: ", energy%core, &
1051 "Hartree energy: ", energy%hartree, &
1052 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1053 END SELECT
1054 ELSE
1055 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1056 "Core Hamiltonian energy: ", energy%core, &
1057 "Hartree energy: ", energy%hartree, &
1058 "Exchange-correlation energy: ", energy%exc + energy%exc_aux_fit
1059 END IF
1060 ELSE
1061!ZMP to print some variables at each step
1062 IF (dft_control%apply_external_density) THEN
1063 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1064 "DOING ZMP CALCULATION FROM EXTERNAL DENSITY "
1065 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1066 "Core Hamiltonian energy: ", energy%core, &
1067 "Hartree energy: ", energy%hartree
1068 ELSE IF (dft_control%apply_external_vxc) THEN
1069 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1070 "DOING ZMP READING EXTERNAL VXC "
1071 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1072 "Core Hamiltonian energy: ", energy%core, &
1073 "Hartree energy: ", energy%hartree
1074 ELSE
1075 IF (psolver .EQ. pw_poisson_implicit) THEN
1076 implicit_ps_ehartree = pw_env%poisson_env%implicit_env%ehartree
1077 bc = pw_env%poisson_env%parameters%ps_implicit_params%boundary_condition
1078 SELECT CASE (bc)
1080 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1081 "Core Hamiltonian energy: ", energy%core, &
1082 "Hartree energy: ", implicit_ps_ehartree, &
1083 "Electric enthalpy: ", energy%hartree, &
1084 "Exchange-correlation energy: ", energy%exc
1085 CASE (periodic_bc, neumann_bc)
1086 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1087 "Core Hamiltonian energy: ", energy%core, &
1088 "Hartree energy: ", energy%hartree, &
1089 "Exchange-correlation energy: ", energy%exc
1090 END SELECT
1091 ELSE
1092 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1093 "Core Hamiltonian energy: ", energy%core, &
1094 "Hartree energy: ", energy%hartree, &
1095 "Exchange-correlation energy: ", energy%exc
1096 END IF
1097 END IF
1098 END IF
1099
1100 IF (dft_control%apply_external_density) THEN
1101 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1102 "Integral of the (density * v_xc): ", energy%exc
1103 END IF
1104
1105 IF (energy%e_hartree /= 0.0_dp) &
1106 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1107 "Coulomb (electron-electron) energy: ", energy%e_hartree
1108 IF (energy%dispersion /= 0.0_dp) &
1109 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1110 "Dispersion energy: ", energy%dispersion
1111 IF (energy%efield /= 0.0_dp) &
1112 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1113 "Electric field interaction energy: ", energy%efield
1114 IF (energy%gcp /= 0.0_dp) &
1115 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1116 "gCP energy: ", energy%gcp
1117 IF (dft_control%qs_control%gapw) THEN
1118 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1119 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit, &
1120 "GAPW| local Eh = 1 center integrals: ", energy%hartree_1c
1121 END IF
1122 IF (dft_control%qs_control%gapw_xc) THEN
1123 WRITE (unit=output_unit, fmt="(/,(T3,A,T61,F20.10))") &
1124 "GAPW| Exc from hard and soft atomic rho1: ", energy%exc1 + energy%exc1_aux_fit
1125 END IF
1126 IF (dft_control%dft_plus_u) THEN
1127 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1128 "DFT+U energy:", energy%dft_plus_u
1129 END IF
1130 IF (qs_env%qmmm) THEN
1131 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1132 "QM/MM Electrostatic energy: ", energy%qmmm_el
1133 IF (qs_env%qmmm_env_qm%image_charge) THEN
1134 WRITE (unit=output_unit, fmt="(T3,A,T61,F20.10)") &
1135 "QM/MM image charge energy: ", energy%image_charge
1136 END IF
1137 END IF
1138 IF (dft_control%qs_control%mulliken_restraint) THEN
1139 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1140 "Mulliken restraint (order_p,energy) : ", mulliken_order_p, energy%mulliken
1141 END IF
1142 IF (dft_control%qs_control%ddapc_restraint) THEN
1143 DO n = 1, SIZE(dft_control%qs_control%ddapc_restraint_control)
1144 ddapc_order_p = &
1145 dft_control%qs_control%ddapc_restraint_control(n)%ddapc_order_p
1146 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1147 "DDAPC restraint (order_p,energy) : ", ddapc_order_p, energy%ddapc_restraint(n)
1148 END DO
1149 END IF
1150 IF (dft_control%qs_control%s2_restraint) THEN
1151 s2_order_p = dft_control%qs_control%s2_restraint_control%s2_order_p
1152 WRITE (unit=output_unit, fmt="(T3,A,T41,2F20.10)") &
1153 "S2 restraint (order_p,energy) : ", s2_order_p, energy%s2_restraint
1154 END IF
1155
1156 END IF ! output_unit
1157 CALL cp_print_key_finished_output(output_unit, logger, input, &
1158 "DFT%SCF%PRINT%DETAILED_ENERGY")
1159
1160 END SUBROUTINE print_detailed_energy
1161
1162! **************************************************************************************************
1163!> \brief compute matrix_vxc, defined via the potential created by qs_vxc_create
1164!> ignores things like tau functional, gapw, sic, ...
1165!> so only OK for GGA & GPW right now
1166!> \param qs_env ...
1167!> \param v_rspace ...
1168!> \param matrix_vxc ...
1169!> \par History
1170!> created 23.10.2012 [Joost VandeVondele]
1171!> \author
1172! **************************************************************************************************
1173 SUBROUTINE compute_matrix_vxc(qs_env, v_rspace, matrix_vxc)
1174 TYPE(qs_environment_type), POINTER :: qs_env
1175 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
1176 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
1177
1178 CHARACTER(LEN=*), PARAMETER :: routinen = 'compute_matrix_vxc'
1179
1180 INTEGER :: handle, ispin
1181 LOGICAL :: gapw
1182 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
1183 TYPE(dft_control_type), POINTER :: dft_control
1184
1185 CALL timeset(routinen, handle)
1186
1187 ! create the matrix using matrix_ks as a template
1188 IF (ASSOCIATED(matrix_vxc)) THEN
1189 CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1190 END IF
1191 CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
1192 ALLOCATE (matrix_vxc(SIZE(matrix_ks)))
1193 DO ispin = 1, SIZE(matrix_ks)
1194 NULLIFY (matrix_vxc(ispin)%matrix)
1195 CALL dbcsr_init_p(matrix_vxc(ispin)%matrix)
1196 CALL dbcsr_copy(matrix_vxc(ispin)%matrix, matrix_ks(ispin)%matrix, &
1197 name="Matrix VXC of spin "//cp_to_string(ispin))
1198 CALL dbcsr_set(matrix_vxc(ispin)%matrix, 0.0_dp)
1199 END DO
1200
1201 ! and integrate
1202 CALL get_qs_env(qs_env, dft_control=dft_control)
1203 gapw = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1204 DO ispin = 1, SIZE(matrix_ks)
1205 CALL integrate_v_rspace(v_rspace=v_rspace(ispin), &
1206 hmat=matrix_vxc(ispin), &
1207 qs_env=qs_env, &
1208 calculate_forces=.false., &
1209 gapw=gapw)
1210 ! scale by the volume element... should really become part of integrate_v_rspace
1211 CALL dbcsr_scale(matrix_vxc(ispin)%matrix, v_rspace(ispin)%pw_grid%dvol)
1212 END DO
1213
1214 CALL timestop(handle)
1215
1216 END SUBROUTINE compute_matrix_vxc
1217
1218! **************************************************************************************************
1219!> \brief Sum up all potentials defined on the grid and integrate
1220!>
1221!> \param qs_env ...
1222!> \param ks_matrix ...
1223!> \param rho ...
1224!> \param my_rho ...
1225!> \param vppl_rspace ...
1226!> \param v_rspace_new ...
1227!> \param v_rspace_new_aux_fit ...
1228!> \param v_tau_rspace ...
1229!> \param v_tau_rspace_aux_fit ...
1230!> \param v_sic_rspace ...
1231!> \param v_spin_ddapc_rest_r ...
1232!> \param v_sccs_rspace ...
1233!> \param v_rspace_embed ...
1234!> \param cdft_control ...
1235!> \param calculate_forces ...
1236!> \par History
1237!> - refactoring 04.03.2011 [MI]
1238!> - SCCS implementation (16.10.2013,MK)
1239!> \author
1240! **************************************************************************************************
1241 SUBROUTINE sum_up_and_integrate(qs_env, ks_matrix, rho, my_rho, &
1242 vppl_rspace, v_rspace_new, &
1243 v_rspace_new_aux_fit, v_tau_rspace, &
1244 v_tau_rspace_aux_fit, &
1245 v_sic_rspace, v_spin_ddapc_rest_r, &
1246 v_sccs_rspace, v_rspace_embed, cdft_control, &
1247 calculate_forces)
1248
1249 TYPE(qs_environment_type), POINTER :: qs_env
1250 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ks_matrix
1251 TYPE(qs_rho_type), POINTER :: rho
1252 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: my_rho
1253 TYPE(pw_r3d_rs_type), POINTER :: vppl_rspace
1254 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new, v_rspace_new_aux_fit, &
1255 v_tau_rspace, v_tau_rspace_aux_fit
1256 TYPE(pw_r3d_rs_type), POINTER :: v_sic_rspace, v_spin_ddapc_rest_r, &
1257 v_sccs_rspace
1258 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_embed
1259 TYPE(cdft_control_type), POINTER :: cdft_control
1260 LOGICAL, INTENT(in) :: calculate_forces
1261
1262 CHARACTER(LEN=*), PARAMETER :: routinen = 'sum_up_and_integrate'
1263
1264 CHARACTER(LEN=default_string_length) :: basis_type
1265 INTEGER :: handle, igroup, ikind, img, ispin, &
1266 nkind, nspins
1267 INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
1268 LOGICAL :: do_ppl, gapw, gapw_xc, lrigpw, rigpw
1269 REAL(kind=dp) :: csign, dvol, fadm
1270 TYPE(admm_type), POINTER :: admm_env
1271 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1272 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ksmat, rho_ao, rho_ao_nokp, smat
1273 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks_aux_fit, &
1274 matrix_ks_aux_fit_dft, rho_ao_aux, &
1275 rho_ao_kp
1276 TYPE(dft_control_type), POINTER :: dft_control
1277 TYPE(kpoint_type), POINTER :: kpoints
1278 TYPE(lri_density_type), POINTER :: lri_density
1279 TYPE(lri_environment_type), POINTER :: lri_env
1280 TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
1281 TYPE(mp_para_env_type), POINTER :: para_env
1282 TYPE(pw_env_type), POINTER :: pw_env
1283 TYPE(pw_poisson_type), POINTER :: poisson_env
1284 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1285 TYPE(pw_r3d_rs_type), POINTER :: v_rspace, vee
1286 TYPE(qs_ks_env_type), POINTER :: ks_env
1287 TYPE(qs_rho_type), POINTER :: rho_aux_fit
1288 TYPE(task_list_type), POINTER :: task_list
1289
1290 CALL timeset(routinen, handle)
1291
1292 NULLIFY (auxbas_pw_pool, dft_control, pw_env, matrix_ks_aux_fit, &
1293 v_rspace, rho_aux_fit, vee, rho_ao, rho_ao_kp, rho_ao_aux, &
1294 ksmat, matrix_ks_aux_fit_dft, lri_env, lri_density, atomic_kind_set, &
1295 rho_ao_nokp, ks_env, admm_env, task_list)
1296
1297 CALL get_qs_env(qs_env, &
1298 dft_control=dft_control, &
1299 pw_env=pw_env, &
1300 v_hartree_rspace=v_rspace, &
1301 vee=vee)
1302
1303 CALL qs_rho_get(rho, rho_ao_kp=rho_ao_kp)
1304 CALL pw_env_get(pw_env, poisson_env=poisson_env, auxbas_pw_pool=auxbas_pw_pool)
1305 gapw = dft_control%qs_control%gapw
1306 gapw_xc = dft_control%qs_control%gapw_xc
1307 do_ppl = dft_control%qs_control%do_ppl_method == do_ppl_grid
1308
1309 rigpw = dft_control%qs_control%rigpw
1310 lrigpw = dft_control%qs_control%lrigpw
1311 IF (lrigpw .OR. rigpw) THEN
1312 CALL get_qs_env(qs_env, &
1313 lri_env=lri_env, &
1314 lri_density=lri_density, &
1315 atomic_kind_set=atomic_kind_set)
1316 END IF
1317
1318 nspins = dft_control%nspins
1319
1320 ! sum up potentials and integrate
1321 IF (ASSOCIATED(v_rspace_new)) THEN
1322 DO ispin = 1, nspins
1323 IF (gapw_xc) THEN
1324 ! SIC not implemented (or at least not tested)
1325 cpassert(dft_control%sic_method_id == sic_none)
1326 !Only the xc potential, because it has to be integrated with the soft basis
1327 CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
1328
1329 ! add the xc part due to v_rspace soft
1330 rho_ao => rho_ao_kp(ispin, :)
1331 ksmat => ks_matrix(ispin, :)
1332 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1333 pmat_kp=rho_ao, hmat_kp=ksmat, &
1334 qs_env=qs_env, &
1335 calculate_forces=calculate_forces, &
1336 gapw=gapw_xc)
1337
1338 ! Now the Hartree potential to be integrated with the full basis
1339 CALL pw_copy(v_rspace, v_rspace_new(ispin))
1340 ELSE
1341 ! Add v_hartree + v_xc = v_rspace_new
1342 CALL pw_axpy(v_rspace, v_rspace_new(ispin), 1.0_dp, v_rspace_new(ispin)%pw_grid%dvol)
1343 END IF ! gapw_xc
1344 IF (dft_control%qs_control%ddapc_explicit_potential) THEN
1345 IF (dft_control%qs_control%ddapc_restraint_is_spin) THEN
1346 IF (ispin == 1) THEN
1347 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1348 ELSE
1349 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), -1.0_dp)
1350 END IF
1351 ELSE
1352 CALL pw_axpy(v_spin_ddapc_rest_r, v_rspace_new(ispin), 1.0_dp)
1353 END IF
1354 END IF
1355 ! CDFT constraint contribution
1356 IF (dft_control%qs_control%cdft) THEN
1357 DO igroup = 1, SIZE(cdft_control%group)
1358 SELECT CASE (cdft_control%group(igroup)%constraint_type)
1360 csign = 1.0_dp
1362 IF (ispin == 1) THEN
1363 csign = 1.0_dp
1364 ELSE
1365 csign = -1.0_dp
1366 END IF
1368 csign = 1.0_dp
1369 IF (ispin == 2) cycle
1371 csign = 1.0_dp
1372 IF (ispin == 1) cycle
1373 CASE DEFAULT
1374 cpabort("Unknown constraint type.")
1375 END SELECT
1376 CALL pw_axpy(cdft_control%group(igroup)%weight, v_rspace_new(ispin), &
1377 csign*cdft_control%strength(igroup))
1378 END DO
1379 END IF
1380 ! functional derivative of the Hartree energy wrt the density in the presence of dielectric
1381 ! (vhartree + v_eps); v_eps is nonzero only if the dielectric constant is defind as a function
1382 ! of the charge density
1383 IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
1384 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1385 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace_new(ispin), dvol)
1386 END IF
1387 ! Add SCCS contribution
1388 IF (dft_control%do_sccs) THEN
1389 CALL pw_axpy(v_sccs_rspace, v_rspace_new(ispin))
1390 END IF
1391 ! External electrostatic potential
1392 IF (dft_control%apply_external_potential) THEN
1393 CALL qmmm_modify_hartree_pot(v_hartree=v_rspace_new(ispin), &
1394 v_qmmm=vee, scale=-1.0_dp)
1395 END IF
1396 IF (do_ppl) THEN
1397 cpassert(.NOT. gapw)
1398 CALL pw_axpy(vppl_rspace, v_rspace_new(ispin), vppl_rspace%pw_grid%dvol)
1399 END IF
1400 ! the electrostatic sic contribution
1401 SELECT CASE (dft_control%sic_method_id)
1402 CASE (sic_none)
1403 !
1405 IF (ispin == 1) THEN
1406 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1407 ELSE
1408 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), 1.0_dp)
1409 END IF
1410 CASE (sic_ad)
1411 CALL pw_axpy(v_sic_rspace, v_rspace_new(ispin), -1.0_dp)
1412 CASE (sic_eo)
1413 ! NOTHING TO BE DONE
1414 END SELECT
1415 ! DFT embedding
1416 IF (dft_control%apply_embed_pot) THEN
1417 CALL pw_axpy(v_rspace_embed(ispin), v_rspace_new(ispin), v_rspace_embed(ispin)%pw_grid%dvol)
1418 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1419 END IF
1420 IF (lrigpw) THEN
1421 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1422 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1423 DO ikind = 1, nkind
1424 lri_v_int(ikind)%v_int = 0.0_dp
1425 END DO
1426 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1427 lri_v_int, calculate_forces, "LRI_AUX")
1428 DO ikind = 1, nkind
1429 CALL para_env%sum(lri_v_int(ikind)%v_int)
1430 END DO
1431 IF (lri_env%exact_1c_terms) THEN
1432 rho_ao => my_rho(ispin, :)
1433 ksmat => ks_matrix(ispin, :)
1434 CALL integrate_v_rspace_diagonal(v_rspace_new(ispin), ksmat(1)%matrix, &
1435 rho_ao(1)%matrix, qs_env, &
1436 calculate_forces, "ORB")
1437 END IF
1438 IF (lri_env%ppl_ri) THEN
1439 CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1440 END IF
1441 ELSEIF (rigpw) THEN
1442 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1443 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1444 DO ikind = 1, nkind
1445 lri_v_int(ikind)%v_int = 0.0_dp
1446 END DO
1447 CALL integrate_v_rspace_one_center(v_rspace_new(ispin), qs_env, &
1448 lri_v_int, calculate_forces, "RI_HXC")
1449 DO ikind = 1, nkind
1450 CALL para_env%sum(lri_v_int(ikind)%v_int)
1451 END DO
1452 ELSE
1453 rho_ao => my_rho(ispin, :)
1454 ksmat => ks_matrix(ispin, :)
1455 CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
1456 pmat_kp=rho_ao, hmat_kp=ksmat, &
1457 qs_env=qs_env, &
1458 calculate_forces=calculate_forces, &
1459 gapw=gapw)
1460 END IF
1461 CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
1462 END DO ! ispin
1463
1464 SELECT CASE (dft_control%sic_method_id)
1465 CASE (sic_none)
1467 CALL auxbas_pw_pool%give_back_pw(v_sic_rspace)
1468 DEALLOCATE (v_sic_rspace)
1469 END SELECT
1470 DEALLOCATE (v_rspace_new)
1471
1472 ELSE
1473 ! not implemented (or at least not tested)
1474 cpassert(dft_control%sic_method_id == sic_none)
1475 cpassert(.NOT. dft_control%qs_control%ddapc_restraint_is_spin)
1476 DO ispin = 1, nspins
1477 ! extra contribution attributed to the dependency of the dielectric constant to the charge density
1478 IF (poisson_env%parameters%solver .EQ. pw_poisson_implicit) THEN
1479 dvol = poisson_env%implicit_env%v_eps%pw_grid%dvol
1480 CALL pw_axpy(poisson_env%implicit_env%v_eps, v_rspace, dvol)
1481 END IF
1482 ! Add SCCS contribution
1483 IF (dft_control%do_sccs) THEN
1484 CALL pw_axpy(v_sccs_rspace, v_rspace)
1485 END IF
1486 ! DFT embedding
1487 IF (dft_control%apply_embed_pot) THEN
1488 CALL pw_axpy(v_rspace_embed(ispin), v_rspace, v_rspace_embed(ispin)%pw_grid%dvol)
1489 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1490 END IF
1491 IF (lrigpw) THEN
1492 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1493 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1494 DO ikind = 1, nkind
1495 lri_v_int(ikind)%v_int = 0.0_dp
1496 END DO
1497 CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
1498 lri_v_int, calculate_forces, "LRI_AUX")
1499 DO ikind = 1, nkind
1500 CALL para_env%sum(lri_v_int(ikind)%v_int)
1501 END DO
1502 IF (lri_env%exact_1c_terms) THEN
1503 rho_ao => my_rho(ispin, :)
1504 ksmat => ks_matrix(ispin, :)
1505 CALL integrate_v_rspace_diagonal(v_rspace, ksmat(1)%matrix, &
1506 rho_ao(1)%matrix, qs_env, &
1507 calculate_forces, "ORB")
1508 END IF
1509 IF (lri_env%ppl_ri) THEN
1510 CALL v_int_ppl_update(qs_env, lri_v_int, calculate_forces)
1511 END IF
1512 ELSEIF (rigpw) THEN
1513 lri_v_int => lri_density%lri_coefs(ispin)%lri_kinds
1514 CALL get_qs_env(qs_env, nkind=nkind, para_env=para_env)
1515 DO ikind = 1, nkind
1516 lri_v_int(ikind)%v_int = 0.0_dp
1517 END DO
1518 CALL integrate_v_rspace_one_center(v_rspace, qs_env, &
1519 lri_v_int, calculate_forces, "RI_HXC")
1520 DO ikind = 1, nkind
1521 CALL para_env%sum(lri_v_int(ikind)%v_int)
1522 END DO
1523 ELSE
1524 rho_ao => my_rho(ispin, :)
1525 ksmat => ks_matrix(ispin, :)
1526 CALL integrate_v_rspace(v_rspace=v_rspace, &
1527 pmat_kp=rho_ao, &
1528 hmat_kp=ksmat, &
1529 qs_env=qs_env, &
1530 calculate_forces=calculate_forces, &
1531 gapw=gapw)
1532 END IF
1533 END DO
1534 END IF ! ASSOCIATED(v_rspace_new)
1535
1536 ! **** LRIGPW: KS matrix from integrated potential
1537 IF (lrigpw) THEN
1538 CALL get_qs_env(qs_env, ks_env=ks_env)
1539 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
1540 CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
1541 DO ispin = 1, nspins
1542 ksmat => ks_matrix(ispin, :)
1543 CALL calculate_lri_ks_matrix(lri_env, lri_v_int, ksmat, atomic_kind_set, &
1544 cell_to_index=cell_to_index)
1545 END DO
1546 IF (calculate_forces) THEN
1547 CALL calculate_lri_forces(lri_env, lri_density, qs_env, rho_ao_kp, atomic_kind_set)
1548 END IF
1549 ELSEIF (rigpw) THEN
1550 CALL get_qs_env(qs_env, matrix_s=smat)
1551 DO ispin = 1, nspins
1552 CALL calculate_ri_ks_matrix(lri_env, lri_v_int, ks_matrix(ispin, 1)%matrix, &
1553 smat(1)%matrix, atomic_kind_set, ispin)
1554 END DO
1555 IF (calculate_forces) THEN
1556 rho_ao_nokp => rho_ao_kp(:, 1)
1557 CALL calculate_ri_forces(lri_env, lri_density, qs_env, rho_ao_nokp, atomic_kind_set)
1558 END IF
1559 END IF
1560
1561 IF (ASSOCIATED(v_tau_rspace)) THEN
1562 IF (lrigpw .OR. rigpw) THEN
1563 cpabort("LRIGPW/RIGPW not implemented for meta-GGAs")
1564 END IF
1565 DO ispin = 1, nspins
1566 CALL pw_scale(v_tau_rspace(ispin), v_tau_rspace(ispin)%pw_grid%dvol)
1567
1568 rho_ao => rho_ao_kp(ispin, :)
1569 ksmat => ks_matrix(ispin, :)
1570 CALL integrate_v_rspace(v_rspace=v_tau_rspace(ispin), &
1571 pmat_kp=rho_ao, hmat_kp=ksmat, &
1572 qs_env=qs_env, &
1573 calculate_forces=calculate_forces, compute_tau=.true., &
1574 gapw=gapw)
1575 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
1576 END DO
1577 DEALLOCATE (v_tau_rspace)
1578 END IF
1579
1580 ! Add contributions from ADMM if requested
1581 IF (dft_control%do_admm) THEN
1582 CALL get_qs_env(qs_env, admm_env=admm_env)
1583 CALL get_admm_env(admm_env, matrix_ks_aux_fit_kp=matrix_ks_aux_fit, rho_aux_fit=rho_aux_fit, &
1584 matrix_ks_aux_fit_dft_kp=matrix_ks_aux_fit_dft)
1585 CALL qs_rho_get(rho_aux_fit, rho_ao_kp=rho_ao_aux)
1586 IF (ASSOCIATED(v_rspace_new_aux_fit)) THEN
1587 DO ispin = 1, nspins
1588 ! Calculate the xc potential
1589 CALL pw_scale(v_rspace_new_aux_fit(ispin), v_rspace_new_aux_fit(ispin)%pw_grid%dvol)
1590
1591 ! set matrix_ks_aux_fit_dft = matrix_ks_aux_fit(k_HF)
1592 DO img = 1, dft_control%nimages
1593 CALL dbcsr_copy(matrix_ks_aux_fit_dft(ispin, img)%matrix, matrix_ks_aux_fit(ispin, img)%matrix, &
1594 name="DFT exch. part of matrix_ks_aux_fit")
1595 END DO
1596
1597 ! Add potential to ks_matrix aux_fit, skip integration if no DFT correction
1598
1599 IF (admm_env%aux_exch_func .NE. do_admm_aux_exch_func_none) THEN
1600
1601 !GPW by default. IF GAPW, then take relevant task list and basis
1602 CALL get_admm_env(admm_env, task_list_aux_fit=task_list)
1603 basis_type = "AUX_FIT"
1604 IF (admm_env%do_gapw) THEN
1605 task_list => admm_env%admm_gapw_env%task_list
1606 basis_type = "AUX_FIT_SOFT"
1607 END IF
1608 fadm = 1.0_dp
1609 ! Calculate bare scaling of force according to Merlot, 1. IF: ADMMP, 2. IF: ADMMS,
1610 IF (admm_env%do_admmp) THEN
1611 fadm = admm_env%gsi(ispin)**2
1612 ELSE IF (admm_env%do_admms) THEN
1613 fadm = (admm_env%gsi(ispin))**(2.0_dp/3.0_dp)
1614 END IF
1615
1616 rho_ao => rho_ao_aux(ispin, :)
1617 ksmat => matrix_ks_aux_fit(ispin, :)
1618
1619 CALL integrate_v_rspace(v_rspace=v_rspace_new_aux_fit(ispin), &
1620 pmat_kp=rho_ao, &
1621 hmat_kp=ksmat, &
1622 qs_env=qs_env, &
1623 calculate_forces=calculate_forces, &
1624 force_adm=fadm, &
1625 gapw=.false., & !even if actual GAPW calculation, want to use AUX_FIT_SOFT
1626 basis_type=basis_type, &
1627 task_list_external=task_list)
1628 END IF
1629
1630 ! matrix_ks_aux_fit_dft(x_DFT)=matrix_ks_aux_fit_dft(old,k_HF)-matrix_ks_aux_fit(k_HF-x_DFT)
1631 DO img = 1, dft_control%nimages
1632 CALL dbcsr_add(matrix_ks_aux_fit_dft(ispin, img)%matrix, &
1633 matrix_ks_aux_fit(ispin, img)%matrix, 1.0_dp, -1.0_dp)
1634 END DO
1635
1636 CALL auxbas_pw_pool%give_back_pw(v_rspace_new_aux_fit(ispin))
1637 END DO
1638 DEALLOCATE (v_rspace_new_aux_fit)
1639 END IF
1640 ! Clean up v_tau_rspace_aux_fit, which is actually not needed
1641 IF (ASSOCIATED(v_tau_rspace_aux_fit)) THEN
1642 DO ispin = 1, nspins
1643 CALL auxbas_pw_pool%give_back_pw(v_tau_rspace_aux_fit(ispin))
1644 END DO
1645 DEALLOCATE (v_tau_rspace_aux_fit)
1646 END IF
1647 END IF
1648
1649 IF (dft_control%apply_embed_pot) DEALLOCATE (v_rspace_embed)
1650
1651 CALL timestop(handle)
1652
1653 END SUBROUTINE sum_up_and_integrate
1654
1655!**************************************************************************
1656!> \brief Calculate the ZMP potential and energy as in Zhao, Morrison Parr
1657!> PRA 50i, 2138 (1994)
1658!> V_c^\lambda defined as int_rho-rho_0/r-r' or rho-rho_0 times a Lagrange
1659!> multiplier, plus Fermi-Amaldi potential that should give the V_xc in the
1660!> limit \lambda --> \infty
1661!>
1662!> \param qs_env ...
1663!> \param v_rspace_new ...
1664!> \param rho ...
1665!> \param exc ...
1666!> \author D. Varsano [daniele.varsano@nano.cnr.it]
1667! **************************************************************************************************
1668 SUBROUTINE calculate_zmp_potential(qs_env, v_rspace_new, rho, exc)
1669
1670 TYPE(qs_environment_type), POINTER :: qs_env
1671 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_new
1672 TYPE(qs_rho_type), POINTER :: rho
1673 REAL(kind=dp) :: exc
1674
1675 CHARACTER(*), PARAMETER :: routinen = 'calculate_zmp_potential'
1676
1677 INTEGER :: handle, my_val, nelectron, nspins
1678 INTEGER, DIMENSION(2) :: nelectron_spin
1679 LOGICAL :: do_zmp_read, fermi_amaldi
1680 REAL(kind=dp) :: lambda
1681 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_ext_r
1682 TYPE(dft_control_type), POINTER :: dft_control
1683 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_ext_g, rho_g
1684 TYPE(pw_env_type), POINTER :: pw_env
1685 TYPE(pw_poisson_type), POINTER :: poisson_env
1686 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1687 TYPE(pw_r3d_rs_type) :: v_xc_rspace
1688 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1689 TYPE(qs_ks_env_type), POINTER :: ks_env
1690 TYPE(section_vals_type), POINTER :: ext_den_section, input
1691
1692!, v_h_gspace, &
1693
1694 CALL timeset(routinen, handle)
1695 NULLIFY (auxbas_pw_pool)
1696 NULLIFY (pw_env)
1697 NULLIFY (poisson_env)
1698 NULLIFY (v_rspace_new)
1699 NULLIFY (dft_control)
1700 NULLIFY (rho_r, rho_g, tot_rho_ext_r, rho_ext_g)
1701 CALL get_qs_env(qs_env=qs_env, &
1702 pw_env=pw_env, &
1703 ks_env=ks_env, &
1704 rho=rho, &
1705 input=input, &
1706 nelectron_spin=nelectron_spin, &
1707 dft_control=dft_control)
1708 CALL pw_env_get(pw_env=pw_env, &
1709 auxbas_pw_pool=auxbas_pw_pool, &
1710 poisson_env=poisson_env)
1711 CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g)
1712 nspins = 1
1713 ALLOCATE (v_rspace_new(nspins))
1714 CALL auxbas_pw_pool%create_pw(pw=v_rspace_new(1))
1715 CALL auxbas_pw_pool%create_pw(pw=v_xc_rspace)
1716
1717 CALL pw_zero(v_rspace_new(1))
1718 do_zmp_read = dft_control%apply_external_vxc
1719 IF (do_zmp_read) THEN
1720 CALL pw_copy(qs_env%external_vxc, v_rspace_new(1))
1721 exc = accurate_dot_product(v_rspace_new(1)%array, rho_r(1)%array)* &
1722 v_rspace_new(1)%pw_grid%dvol
1723 ELSE
1724 block
1725 REAL(kind=dp) :: factor
1726 TYPE(pw_c1d_gs_type) :: rho_eff_gspace, v_xc_gspace
1727 CALL auxbas_pw_pool%create_pw(pw=rho_eff_gspace)
1728 CALL auxbas_pw_pool%create_pw(pw=v_xc_gspace)
1729 CALL pw_zero(rho_eff_gspace)
1730 CALL pw_zero(v_xc_gspace)
1731 CALL pw_zero(v_xc_rspace)
1732 factor = pw_integrate_function(rho_g(1))
1733 CALL qs_rho_get(qs_env%rho_external, &
1734 rho_g=rho_ext_g, &
1735 tot_rho_r=tot_rho_ext_r)
1736 factor = tot_rho_ext_r(1)/factor
1737
1738 CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1739 CALL pw_axpy(rho_ext_g(1), rho_eff_gspace, alpha=-1.0_dp)
1740 ext_den_section => section_vals_get_subs_vals(input, "DFT%EXTERNAL_DENSITY")
1741 CALL section_vals_val_get(ext_den_section, "LAMBDA", r_val=lambda)
1742 CALL section_vals_val_get(ext_den_section, "ZMP_CONSTRAINT", i_val=my_val)
1743 CALL section_vals_val_get(ext_den_section, "FERMI_AMALDI", l_val=fermi_amaldi)
1744
1745 CALL pw_scale(rho_eff_gspace, a=lambda)
1746 nelectron = nelectron_spin(1)
1747 factor = -1.0_dp/nelectron
1748 CALL pw_axpy(rho_g(1), rho_eff_gspace, alpha=factor)
1749
1750 CALL pw_poisson_solve(poisson_env, rho_eff_gspace, vhartree=v_xc_gspace)
1751 CALL pw_transfer(v_xc_gspace, v_rspace_new(1))
1752 CALL pw_copy(v_rspace_new(1), v_xc_rspace)
1753
1754 exc = 0.0_dp
1755 exc = pw_integral_ab(v_rspace_new(1), rho_r(1))
1756
1757!Note that this is not the xc energy but \int(\rho*v_xc)
1758!Vxc---> v_rspace_new
1759!Exc---> energy%exc
1760 CALL auxbas_pw_pool%give_back_pw(rho_eff_gspace)
1761 CALL auxbas_pw_pool%give_back_pw(v_xc_gspace)
1762 END block
1763 END IF
1764
1765 CALL auxbas_pw_pool%give_back_pw(v_xc_rspace)
1766
1767 CALL timestop(handle)
1768
1769 END SUBROUTINE calculate_zmp_potential
1770
1771! **************************************************************************************************
1772!> \brief ...
1773!> \param qs_env ...
1774!> \param rho ...
1775!> \param v_rspace_embed ...
1776!> \param dft_control ...
1777!> \param embed_corr ...
1778!> \param just_energy ...
1779! **************************************************************************************************
1780 SUBROUTINE get_embed_potential_energy(qs_env, rho, v_rspace_embed, dft_control, embed_corr, &
1781 just_energy)
1782 TYPE(qs_environment_type), POINTER :: qs_env
1783 TYPE(qs_rho_type), POINTER :: rho
1784 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_rspace_embed
1785 TYPE(dft_control_type), POINTER :: dft_control
1786 REAL(kind=dp) :: embed_corr
1787 LOGICAL :: just_energy
1788
1789 CHARACTER(*), PARAMETER :: routinen = 'get_embed_potential_energy'
1790
1791 INTEGER :: handle, ispin
1792 REAL(kind=dp) :: embed_corr_local
1793 TYPE(pw_env_type), POINTER :: pw_env
1794 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1795 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
1796
1797 CALL timeset(routinen, handle)
1798
1799 NULLIFY (auxbas_pw_pool)
1800 NULLIFY (pw_env)
1801 NULLIFY (rho_r)
1802 CALL get_qs_env(qs_env=qs_env, &
1803 pw_env=pw_env, &
1804 rho=rho)
1805 CALL pw_env_get(pw_env=pw_env, &
1806 auxbas_pw_pool=auxbas_pw_pool)
1807 CALL qs_rho_get(rho, rho_r=rho_r)
1808 ALLOCATE (v_rspace_embed(dft_control%nspins))
1809
1810 embed_corr = 0.0_dp
1811
1812 DO ispin = 1, dft_control%nspins
1813 CALL auxbas_pw_pool%create_pw(pw=v_rspace_embed(ispin))
1814 CALL pw_zero(v_rspace_embed(ispin))
1815
1816 CALL pw_copy(qs_env%embed_pot, v_rspace_embed(ispin))
1817 embed_corr_local = 0.0_dp
1818
1819 ! Spin embedding potential in open-shell case
1820 IF (dft_control%nspins .EQ. 2) THEN
1821 IF (ispin .EQ. 1) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), 1.0_dp)
1822 IF (ispin .EQ. 2) CALL pw_axpy(qs_env%spin_embed_pot, v_rspace_embed(ispin), -1.0_dp)
1823 END IF
1824 ! Integrate the density*potential
1825 embed_corr_local = pw_integral_ab(v_rspace_embed(ispin), rho_r(ispin))
1826
1827 embed_corr = embed_corr + embed_corr_local
1828
1829 END DO
1830
1831 ! If only energy requiested we delete the potential
1832 IF (just_energy) THEN
1833 DO ispin = 1, dft_control%nspins
1834 CALL auxbas_pw_pool%give_back_pw(v_rspace_embed(ispin))
1835 END DO
1836 DEALLOCATE (v_rspace_embed)
1837 END IF
1838
1839 CALL timestop(handle)
1840
1841 END SUBROUTINE get_embed_potential_energy
1842
1843END 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...
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)
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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
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)
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_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, 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:1449
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:509
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.