(git:9111030)
Loading...
Searching...
No Matches
xc_pot_saop.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Calculate the saop potential
10! **************************************************************************************************
11MODULE xc_pot_saop
12 USE atomic_kind_types, ONLY: atomic_kind_type,&
13 get_atomic_kind
14 USE basis_set_types, ONLY: gto_basis_set_type
15 USE cp_array_utils, ONLY: cp_1d_r_p_type
16 USE cp_control_types, ONLY: dft_control_type
17 USE cp_dbcsr_api, ONLY: dbcsr_copy,&
18 dbcsr_deallocate_matrix,&
19 dbcsr_p_type,&
20 dbcsr_set
21 USE cp_dbcsr_operations, ONLY: cp_dbcsr_plus_fm_fm_t,&
22 dbcsr_allocate_matrix_set,&
23 dbcsr_deallocate_matrix_set
24 USE cp_fm_types, ONLY: cp_fm_create,&
25 cp_fm_get_info,&
26 cp_fm_get_submatrix,&
27 cp_fm_p_type,&
28 cp_fm_release,&
29 cp_fm_set_all,&
30 cp_fm_set_submatrix,&
31 cp_fm_type
32 USE input_constants, ONLY: do_method_gapw,&
33 oe_gllb,&
34 oe_lb,&
35 oe_saop,&
36 xc_funct_no_shortcut
37 USE input_section_types, ONLY: &
38 section_vals_create, section_vals_duplicate, section_vals_get_subs_vals, &
39 section_vals_release, section_vals_retain, section_vals_set_subs_vals, section_vals_type, &
40 section_vals_val_get, section_vals_val_set
41 USE kinds, ONLY: dp
42 USE mathconstants, ONLY: pi
43 USE message_passing, ONLY: mp_para_env_type
44 USE pw_env_types, ONLY: pw_env_get,&
45 pw_env_type
46 USE pw_methods, ONLY: pw_axpy,&
47 pw_copy,&
48 pw_scale,&
49 pw_zero
50 USE pw_pool_types, ONLY: pw_pool_type
51 USE pw_types, ONLY: pw_c1d_gs_type,&
52 pw_r3d_rs_type
53 USE qs_collocate_density, ONLY: calculate_rho_elec
54 USE qs_environment_types, ONLY: get_qs_env,&
55 qs_environment_type
56 USE qs_gapw_densities, ONLY: prepare_gapw_den
57 USE qs_grid_atom, ONLY: grid_atom_type
58 USE qs_harmonics_atom, ONLY: harmonics_atom_type
59 USE qs_integrate_potential, ONLY: integrate_v_rspace
60 USE qs_kind_types, ONLY: get_qs_kind,&
61 qs_kind_type
62 USE qs_ks_atom, ONLY: update_ks_atom
63 USE qs_ks_types, ONLY: qs_ks_env_type
64 USE qs_local_rho_types, ONLY: local_rho_set_create,&
65 local_rho_set_release,&
66 local_rho_type
67 USE qs_mo_types, ONLY: get_mo_set,&
68 mo_set_type
69 USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
70 USE qs_oce_types, ONLY: oce_matrix_type
71 USE qs_rho_atom_methods, ONLY: allocate_rho_atom_internals,&
72 calculate_rho_atom_coeff
73 USE qs_rho_atom_types, ONLY: get_rho_atom,&
74 rho_atom_coeff,&
75 rho_atom_type
76 USE qs_rho_types, ONLY: qs_rho_get,&
77 qs_rho_type
78 USE qs_vxc_atom, ONLY: calc_rho_angular,&
79 gavxcgb_nogc
80 USE util, ONLY: get_limit
81 USE virial_types, ONLY: virial_type
82 USE xc, ONLY: xc_vxc_pw_create
83 USE xc_atom, ONLY: fill_rho_set,&
84 vxc_of_r_new,&
85 xc_rho_set_atom_update
86 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
87 xc_dset_create,&
88 xc_dset_get_derivative,&
89 xc_dset_release,&
90 xc_dset_zero_all
91 USE xc_derivative_types, ONLY: xc_derivative_get,&
92 xc_derivative_type
93 USE xc_derivatives, ONLY: xc_functionals_eval
94 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_setall,&
95 xc_rho_cflags_type
96 USE xc_rho_set_types, ONLY: xc_rho_set_create,&
97 xc_rho_set_release,&
98 xc_rho_set_type,&
99 xc_rho_set_update
100 USE xc_xbecke88, ONLY: xb88_lda_info,&
101 xb88_lsd_info
102#include "./base/base_uses.f90"
103
104 IMPLICIT NONE
105
106 PRIVATE
107
108 PUBLIC :: add_saop_pot
109
110 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pot_saop'
111
112 ! should be eliminated
113 REAL(KIND=dp), PARAMETER :: alpha = 1.19_dp, beta = 0.01_dp, k_rho = 0.42_dp
114 REAL(KIND=dp), PARAMETER :: kappa = 0.804_dp, mu = 0.21951_dp, &
115 beta_ec = 0.066725_dp, gamma_saop = 0.031091_dp
116
117CONTAINS
118
119! **************************************************************************************************
120!> \brief ...
121!> \param ks_matrix ...
122!> \param qs_env ...
123!> \param oe_corr ...
124! **************************************************************************************************
125 SUBROUTINE add_saop_pot(ks_matrix, qs_env, oe_corr)
126
127 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
128 TYPE(qs_environment_type), POINTER :: qs_env
129 INTEGER, INTENT(IN) :: oe_corr
130
131 INTEGER :: dft_method_id, homo, i, ispin, j, k, &
132 nspins, orb, xc_deriv_method_id, &
133 xc_rho_smooth_id
134 INTEGER, DIMENSION(2) :: ncol, nrow
135 INTEGER, DIMENSION(2, 3) :: bo
136 LOGICAL :: compute_virial, gapw, lsd
137 REAL(kind=dp) :: density_cut, efac, gradient_cut, &
138 tau_cut, we_gllb, we_lb, xc_energy
139 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_col
140 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
141 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
142 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: e_uniform
143 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: single_mo_coeff
144 TYPE(cp_fm_type), POINTER :: mo_coeff
145 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: orbital_density_matrix, rho_struct_ao
146 TYPE(mo_set_type), DIMENSION(:), POINTER :: molecular_orbitals
147 TYPE(pw_c1d_gs_type) :: orbital_g
148 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
149 TYPE(pw_env_type), POINTER :: pw_env
150 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
151 TYPE(pw_r3d_rs_type) :: orbital
152 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: vxc_gllb, vxc_saop
153 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_struct_r, tau, vxc_lb, &
154 vxc_tau, vxc_tmp
155 TYPE(pw_r3d_rs_type), POINTER :: weights
156 TYPE(qs_ks_env_type), POINTER :: ks_env
157 TYPE(qs_rho_type), POINTER :: rho_struct
158 TYPE(section_vals_type), POINTER :: input, xc_fun_section_orig, &
159 xc_fun_section_tmp, xc_section_orig, &
160 xc_section_tmp
161 TYPE(virial_type), POINTER :: virial
162 TYPE(xc_derivative_set_type) :: deriv_set
163 TYPE(xc_derivative_type), POINTER :: deriv
164 TYPE(xc_rho_cflags_type) :: needs
165 TYPE(xc_rho_set_type) :: rho_set
166
167 NULLIFY (ks_env, pw_env, auxbas_pw_pool, input)
168 NULLIFY (rho_g, rho_r, tau, rho_struct, e_uniform)
169 NULLIFY (vxc_lb, vxc_tmp, vxc_tau)
170 NULLIFY (mo_eigenvalues, deriv, rho_struct_r, rho_struct_ao)
171 NULLIFY (orbital_density_matrix, xc_section_tmp, xc_fun_section_tmp)
172
173 CALL get_qs_env(qs_env, &
174 ks_env=ks_env, &
175 rho=rho_struct, &
176 xcint_weights=weights, &
177 pw_env=pw_env, &
178 input=input, &
179 virial=virial, &
180 mos=molecular_orbitals)
181 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
182 CALL section_vals_val_get(input, "DFT%QS%METHOD", i_val=dft_method_id)
183 gapw = (dft_method_id == do_method_gapw)
184
185 xc_section_orig => section_vals_get_subs_vals(input, "DFT%XC")
186 CALL section_vals_retain(xc_section_orig)
187 CALL section_vals_duplicate(xc_section_orig, xc_section_tmp)
188
189 CALL section_vals_val_get(xc_section_orig, "DENSITY_CUTOFF", &
190 r_val=density_cut)
191 CALL section_vals_val_get(xc_section_orig, "GRADIENT_CUTOFF", &
192 r_val=gradient_cut)
193 CALL section_vals_val_get(xc_section_orig, "TAU_CUTOFF", &
194 r_val=tau_cut)
195
196 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
197
198 CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
199 IF (lsd) THEN
200 nspins = 2
201 ELSE
202 nspins = 1
203 END IF
204
205 ALLOCATE (single_mo_coeff(nspins))
206 CALL dbcsr_allocate_matrix_set(orbital_density_matrix, nspins)
207 CALL qs_rho_get(rho_struct, rho_r=rho_struct_r, rho_ao=rho_struct_ao)
208 rho_r => rho_struct_r
209 DO ispin = 1, nspins
210 ALLOCATE (orbital_density_matrix(ispin)%matrix)
211 CALL dbcsr_copy(orbital_density_matrix(ispin)%matrix, &
212 rho_struct_ao(ispin)%matrix, "orbital density")
213 END DO
214 bo = rho_r(1)%pw_grid%bounds_local
215
216 !---------------------------!
217 ! create the density needed !
218 !---------------------------!
219 CALL xc_rho_set_create(rho_set, bo, &
220 density_cut, &
221 gradient_cut, &
222 tau_cut)
223 CALL xc_rho_cflags_setall(needs, .false.)
224 IF (lsd) THEN
225 CALL xb88_lsd_info(needs=needs)
226 needs%norm_drho = .true.
227 ELSE
228 CALL xb88_lda_info(needs=needs)
229 END IF
230 CALL section_vals_val_get(xc_section_orig, "XC_GRID%XC_DERIV", &
231 i_val=xc_deriv_method_id)
232 CALL section_vals_val_get(xc_section_orig, "XC_GRID%XC_SMOOTH_RHO", &
233 i_val=xc_rho_smooth_id)
234 CALL xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, &
235 xc_deriv_method_id, &
236 xc_rho_smooth_id, &
237 auxbas_pw_pool)
238
239 !----------------------------------------!
240 ! Construct the LB94 potential in vxc_LB !
241 !----------------------------------------!
242 xc_fun_section_orig => section_vals_get_subs_vals(xc_section_orig, &
243 "XC_FUNCTIONAL")
244 CALL section_vals_create(xc_fun_section_tmp, xc_fun_section_orig%section)
245 CALL section_vals_val_set(xc_fun_section_tmp, "_SECTION_PARAMETERS_", &
246 i_val=xc_funct_no_shortcut)
247 CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
248 l_val=.true.)
249 CALL section_vals_set_subs_vals(xc_section_tmp, "XC_FUNCTIONAL", &
250 xc_fun_section_tmp)
251
252 cpassert(.NOT. compute_virial)
253 CALL xc_vxc_pw_create(vxc_tmp, vxc_tau, xc_energy, rho_r, rho_g, tau, &
254 xc_section_tmp, weights, auxbas_pw_pool, &
255 compute_virial=.false., virial_xc=virial_xc_tmp)
256
257 CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
258 l_val=.false.)
259 CALL section_vals_val_set(xc_fun_section_tmp, "PZ81%_SECTION_PARAMETERS_", &
260 l_val=.true.)
261
262 cpassert(.NOT. compute_virial)
263 CALL xc_vxc_pw_create(vxc_lb, vxc_tau, xc_energy, rho_r, rho_g, tau, &
264 xc_section_tmp, weights, auxbas_pw_pool, &
265 compute_virial=.false., virial_xc=virial_xc_tmp)
266
267 DO ispin = 1, nspins
268 CALL pw_axpy(vxc_tmp(ispin), vxc_lb(ispin), alpha)
269 END DO
270
271 DO ispin = 1, nspins
272 CALL add_lb_pot(vxc_tmp(ispin)%array, rho_set, lsd, ispin)
273 CALL pw_axpy(vxc_tmp(ispin), vxc_lb(ispin), -1.0_dp)
274 END DO
275
276 !-----------------------------------------------------------------------------------!
277 ! Construct 2 times PBE one particle density from the PZ correlation energy density !
278 !-----------------------------------------------------------------------------------!
279 CALL xc_dset_create(deriv_set, local_bounds=bo)
280 CALL xc_functionals_eval(xc_fun_section_tmp, &
281 lsd=lsd, &
282 rho_set=rho_set, &
283 deriv_set=deriv_set, &
284 deriv_order=0)
285
286 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
287 CALL xc_derivative_get(deriv, deriv_data=e_uniform)
288
289 ALLOCATE (vxc_gllb(nspins))
290 DO ispin = 1, nspins
291 CALL auxbas_pw_pool%create_pw(vxc_gllb(ispin))
292 END DO
293
294 DO ispin = 1, nspins
295 CALL calc_2excpbe(vxc_gllb(ispin)%array, rho_set, e_uniform, lsd)
296 END DO
297
298 CALL xc_dset_release(deriv_set)
299
300 CALL auxbas_pw_pool%create_pw(orbital)
301 CALL auxbas_pw_pool%create_pw(orbital_g)
302
303 DO ispin = 1, nspins
304
305 CALL get_mo_set(molecular_orbitals(ispin), &
306 mo_coeff=mo_coeff, &
307 eigenvalues=mo_eigenvalues, &
308 homo=homo)
309 CALL cp_fm_create(single_mo_coeff(ispin), &
310 mo_coeff%matrix_struct, &
311 "orbital density matrix")
312
313 CALL cp_fm_get_info(single_mo_coeff(ispin), &
314 nrow_global=nrow(ispin), ncol_global=ncol(ispin))
315 ALLOCATE (coeff_col(nrow(ispin), 1))
316
317 CALL pw_zero(vxc_tmp(ispin))
318
319 DO orb = 1, homo - 1
320
321 efac = k_rho*sqrt(mo_eigenvalues(homo) - mo_eigenvalues(orb))
322 IF (.NOT. lsd) efac = 2.0_dp*efac
323
324 CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
325 CALL cp_fm_get_submatrix(mo_coeff, coeff_col, &
326 1, orb, nrow(ispin), 1)
327 CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
328 1, orb)
329 CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
330 CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
331 matrix_v=single_mo_coeff(ispin), &
332 ncol=ncol(ispin), &
333 alpha=1.0_dp)
334 CALL pw_zero(orbital)
335 CALL pw_zero(orbital_g)
336 CALL calculate_rho_elec(matrix_p=orbital_density_matrix(ispin)%matrix, &
337 rho=orbital, rho_gspace=orbital_g, &
338 ks_env=ks_env)
339
340 CALL pw_axpy(orbital, vxc_tmp(ispin), efac)
341
342 END DO
343 DEALLOCATE (coeff_col)
344
345 DO k = bo(1, 3), bo(2, 3)
346 DO j = bo(1, 2), bo(2, 2)
347 DO i = bo(1, 1), bo(2, 1)
348 IF (rho_r(ispin)%array(i, j, k) > density_cut) THEN
349 vxc_tmp(ispin)%array(i, j, k) = vxc_tmp(ispin)%array(i, j, k)/ &
350 rho_r(ispin)%array(i, j, k)
351 ELSE
352 vxc_tmp(ispin)%array(i, j, k) = 0.0_dp
353 END IF
354 END DO
355 END DO
356 END DO
357
358 CALL pw_axpy(vxc_tmp(ispin), vxc_gllb(ispin), 1.0_dp)
359
360 END DO
361
362 !---------------!
363 ! Assemble SAOP !
364 !---------------!
365 ALLOCATE (vxc_saop(nspins))
366
367 DO ispin = 1, nspins
368
369 CALL get_mo_set(molecular_orbitals(ispin), &
370 mo_coeff=mo_coeff, &
371 eigenvalues=mo_eigenvalues, &
372 homo=homo)
373 CALL auxbas_pw_pool%create_pw(vxc_saop(ispin))
374 CALL pw_zero(vxc_saop(ispin))
375
376 ALLOCATE (coeff_col(nrow(ispin), 1))
377
378 DO orb = 1, homo
379
380 we_lb = exp(-2.0_dp*(mo_eigenvalues(homo) - mo_eigenvalues(orb))**2)
381 we_gllb = 1.0_dp - we_lb
382 IF (.NOT. lsd) THEN
383 we_lb = 2.0_dp*we_lb
384 we_gllb = 2.0_dp*we_gllb
385 END IF
386
387 vxc_tmp(ispin)%array = we_lb*vxc_lb(ispin)%array + &
388 we_gllb*vxc_gllb(ispin)%array
389
390 CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
391 CALL cp_fm_get_submatrix(mo_coeff, coeff_col, &
392 1, orb, nrow(ispin), 1)
393 CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
394 1, orb)
395 CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
396 CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
397 matrix_v=single_mo_coeff(ispin), &
398 ncol=ncol(ispin), &
399 alpha=1.0_dp)
400 CALL pw_zero(orbital)
401 CALL pw_zero(orbital_g)
402 CALL calculate_rho_elec(matrix_p=orbital_density_matrix(ispin)%matrix, &
403 rho=orbital, rho_gspace=orbital_g, &
404 ks_env=ks_env)
405
406 vxc_saop(ispin)%array = vxc_saop(ispin)%array + &
407 orbital%array*vxc_tmp(ispin)%array
408
409 END DO
410
411 CALL dbcsr_deallocate_matrix(orbital_density_matrix(ispin)%matrix)
412
413 DEALLOCATE (coeff_col)
414
415 DO k = bo(1, 3), bo(2, 3)
416 DO j = bo(1, 2), bo(2, 2)
417 DO i = bo(1, 1), bo(2, 1)
418 IF (rho_r(ispin)%array(i, j, k) > density_cut) THEN
419 vxc_saop(ispin)%array(i, j, k) = vxc_saop(ispin)%array(i, j, k)/ &
420 rho_r(ispin)%array(i, j, k)
421 ELSE
422 vxc_saop(ispin)%array(i, j, k) = 0.0_dp
423 END IF
424 END DO
425 END DO
426 END DO
427
428 END DO
429
430 CALL cp_fm_release(single_mo_coeff)
431
432 CALL xc_rho_set_release(rho_set, auxbas_pw_pool)
433 CALL auxbas_pw_pool%give_back_pw(orbital)
434 CALL auxbas_pw_pool%give_back_pw(orbital_g)
435
436 !--------------------!
437 ! Do the integration !
438 !--------------------!
439 DO ispin = 1, nspins
440
441 IF (oe_corr == oe_lb) THEN
442 CALL pw_copy(vxc_lb(ispin), vxc_saop(ispin))
443 ELSE IF (oe_corr == oe_gllb) THEN
444 CALL pw_copy(vxc_gllb(ispin), vxc_saop(ispin))
445 END IF
446 CALL pw_scale(vxc_saop(ispin), vxc_saop(ispin)%pw_grid%dvol)
447
448 CALL integrate_v_rspace(v_rspace=vxc_saop(ispin), pmat=rho_struct_ao(ispin), &
449 hmat=ks_matrix(ispin), qs_env=qs_env, &
450 calculate_forces=.false., &
451 gapw=gapw)
452
453 END DO
454
455 DO ispin = 1, nspins
456 CALL auxbas_pw_pool%give_back_pw(vxc_saop(ispin))
457 CALL auxbas_pw_pool%give_back_pw(vxc_gllb(ispin))
458 CALL vxc_lb(ispin)%release()
459 CALL vxc_tmp(ispin)%release()
460 END DO
461 DEALLOCATE (vxc_gllb, vxc_lb, vxc_tmp, orbital_density_matrix)
462
463 DEALLOCATE (vxc_saop)
464
465 CALL section_vals_release(xc_fun_section_tmp)
466 CALL section_vals_release(xc_section_tmp)
467 CALL section_vals_release(xc_section_orig)
468
469 !-----------------------!
470 ! Call the GAPW routine !
471 !-----------------------!
472 IF (gapw) THEN
473 CALL gapw_add_atomic_saop_pot(qs_env, oe_corr)
474 END IF
475
476 END SUBROUTINE add_saop_pot
477
478! **************************************************************************************************
479!> \brief ...
480!> \param qs_env ...
481!> \param oe_corr ...
482! **************************************************************************************************
483 SUBROUTINE gapw_add_atomic_saop_pot(qs_env, oe_corr)
484
485 TYPE(qs_environment_type), POINTER :: qs_env
486 INTEGER, INTENT(IN) :: oe_corr
487
488 INTEGER :: ia, iat, iatom, ikind, ir, ispin, na, &
489 natom, nr, ns, nspins, orb
490 INTEGER, DIMENSION(2) :: bo, homo, ncol, nrow
491 INTEGER, DIMENSION(2, 3) :: bounds
492 INTEGER, DIMENSION(:), POINTER :: atom_list
493 LOGICAL :: accint, lsd, paw_atom
494 REAL(dp), DIMENSION(:, :, :), POINTER :: tau
495 REAL(kind=dp) :: agr, alpha, density_cut, efac, exc, &
496 gradient_cut, tau_cut, we_gllb, we_lb
497 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_col, weight_h, weight_s
498 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_uniform, rho_h, rho_s, vtau, &
499 vxc_gllb_h, vxc_gllb_s, vxc_lb_h, vxc_lb_s, vxc_saop_h, vxc_saop_s, vxc_tmp_h, vxc_tmp_s
500 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg
501 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
502 TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: mo_eigenvalues
503 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff
504 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: single_mo_coeff
505 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, orbital_density_matrix, &
506 rho_struct_ao
507 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat
508 TYPE(dft_control_type), POINTER :: dft_control
509 TYPE(grid_atom_type), POINTER :: atomic_grid, grid_atom
510 TYPE(gto_basis_set_type), POINTER :: orb_basis
511 TYPE(harmonics_atom_type), POINTER :: harmonics
512 TYPE(local_rho_type), POINTER :: local_rho_set
513 TYPE(mo_set_type), DIMENSION(:), POINTER :: molecular_orbitals
514 TYPE(mp_para_env_type), POINTER :: para_env
515 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
516 POINTER :: sab
517 TYPE(oce_matrix_type), POINTER :: oce
518 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
519 TYPE(qs_rho_type), POINTER :: rho_structure
520 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
521 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
522 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
523 TYPE(rho_atom_type), POINTER :: rho_atom
524 TYPE(section_vals_type), POINTER :: input, xc_fun_section_orig, &
525 xc_fun_section_tmp, xc_section_orig, &
526 xc_section_tmp
527 TYPE(xc_derivative_set_type) :: deriv_set
528 TYPE(xc_derivative_type), POINTER :: deriv
529 TYPE(xc_rho_cflags_type) :: needs, needs_orbs
530 TYPE(xc_rho_set_type) :: orb_rho_set_h, orb_rho_set_s, rho_set_h, &
531 rho_set_s
532
533 NULLIFY (rho_h, rho_s, vxc_lb_h, vxc_lb_s, vxc_gllb_h, vxc_gllb_s, &
534 vxc_tmp_h, vxc_tmp_s, vtau, dummy, e_uniform, drho_h, drho_s, vxg, atom_list, &
535 atomic_kind_set, qs_kind_set, deriv, atomic_grid, rho_struct_ao, &
536 harmonics, molecular_orbitals, rho_structure, r_h, r_s, dr_h, dr_s, &
537 r_h_d, r_s_d, rho_atom_set, rho_atom, para_env, &
538 mo_eigenvalues, local_rho_set, matrix_ks, &
539 orbital_density_matrix, vxc_saop_h, vxc_saop_s)
540
541 ! tau is needed for fill_rho_set, but should never be used
542 NULLIFY (tau)
543 NULLIFY (dft_control, oce, sab)
544
545 CALL get_qs_env(qs_env, input=input, &
546 rho=rho_structure, &
547 mos=molecular_orbitals, &
548 atomic_kind_set=atomic_kind_set, &
549 qs_kind_set=qs_kind_set, &
550 rho_atom_set=rho_atom_set, &
551 matrix_ks=matrix_ks, &
552 dft_control=dft_control, &
553 para_env=para_env, &
554 oce=oce, sab_orb=sab)
555
556 CALL qs_rho_get(rho_structure, rho_ao=rho_struct_ao)
557
558 xc_section_orig => section_vals_get_subs_vals(input, "DFT%XC")
559 CALL section_vals_retain(xc_section_orig)
560 CALL section_vals_duplicate(xc_section_orig, xc_section_tmp)
561
562 accint = dft_control%qs_control%gapw_control%accurate_xcint
563
564 ! [SC] the following code can be traced back to SVN rev. 4296 (git:f97138b) that
565 ! has removed the component 'nspins' from the derived type 'dft_control_type'.
566 ! Is it worth to remove the code below in favour of 'dft_control%nspins'
567 ! since its reintroduction? Note that in case of ROKS calculations,
568 ! 'lsd == .FALSE.' but 'dft_control%nspins == 2'.
569 CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
570 IF (lsd) THEN
571 nspins = 2
572 ELSE
573 nspins = 1
574 END IF
575
576 CALL section_vals_val_get(xc_section_orig, "DENSITY_CUTOFF", &
577 r_val=density_cut)
578 CALL section_vals_val_get(xc_section_orig, "GRADIENT_CUTOFF", &
579 r_val=gradient_cut)
580 CALL section_vals_val_get(xc_section_orig, "TAU_CUTOFF", &
581 r_val=tau_cut)
582
583 ! remap pointer
584 ns = SIZE(rho_struct_ao)
585 psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
586 CALL calculate_rho_atom_coeff(qs_env, psmat, rho_atom_set, qs_kind_set, oce, sab, para_env)
587 CALL prepare_gapw_den(qs_env)
588
589 ALLOCATE (mo_coeff(nspins), single_mo_coeff(nspins), mo_eigenvalues(nspins))
590
591 CALL dbcsr_allocate_matrix_set(orbital_density_matrix, nspins)
592
593 DO ispin = 1, nspins
594 CALL get_mo_set(molecular_orbitals(ispin), &
595 mo_coeff=mo_coeff(ispin)%matrix, &
596 eigenvalues=mo_eigenvalues(ispin)%array, &
597 homo=homo(ispin))
598 CALL cp_fm_create(single_mo_coeff(ispin), &
599 mo_coeff(ispin)%matrix%matrix_struct, &
600 "orbital density matrix")
601 CALL cp_fm_get_info(single_mo_coeff(ispin), &
602 nrow_global=nrow(ispin), ncol_global=ncol(ispin))
603 ALLOCATE (orbital_density_matrix(ispin)%matrix)
604 CALL dbcsr_copy(orbital_density_matrix(ispin)%matrix, &
605 rho_struct_ao(ispin)%matrix, &
606 "orbital density")
607 END DO
608 CALL local_rho_set_create(local_rho_set)
609 CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
610 qs_kind_set, dft_control, para_env)
611
612 DO ikind = 1, SIZE(atomic_kind_set)
613 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
614
615 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
616 harmonics=harmonics, grid_atom=atomic_grid)
617 IF (.NOT. paw_atom) cycle
618
619 nr = atomic_grid%nr
620 na = atomic_grid%ng_sphere
621 bounds(1:2, 1:3) = 1
622 bounds(2, 1) = na
623 bounds(2, 2) = nr
624
625 CALL xc_dset_create(deriv_set, local_bounds=bounds)
626
627 CALL xc_rho_set_create(rho_set_h, bounds, density_cut, &
628 gradient_cut, tau_cut)
629 CALL xc_rho_set_create(rho_set_s, bounds, density_cut, &
630 gradient_cut, tau_cut)
631 CALL xc_rho_set_create(orb_rho_set_h, bounds, density_cut, &
632 gradient_cut, tau_cut)
633 CALL xc_rho_set_create(orb_rho_set_s, bounds, density_cut, &
634 gradient_cut, tau_cut)
635
636 CALL xc_rho_cflags_setall(needs, .false.)
637 IF (lsd) THEN
638 CALL xb88_lsd_info(needs=needs)
639 needs%norm_drho = .true.
640 ELSE
641 CALL xb88_lda_info(needs=needs)
642 END IF
643 CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
644 CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
645 CALL xc_rho_cflags_setall(needs_orbs, .false.)
646 needs_orbs%rho = .true.
647 IF (lsd) needs_orbs%rho_spin = .true.
648 CALL xc_rho_set_atom_update(orb_rho_set_h, needs, nspins, bounds)
649 CALL xc_rho_set_atom_update(orb_rho_set_s, needs, nspins, bounds)
650
651 ALLOCATE (rho_h(1:na, 1:nr, 1:nspins), rho_s(1:na, 1:nr, 1:nspins))
652 ALLOCATE (weight_h(1:na, 1:nr), weight_s(1:na, 1:nr))
653 ALLOCATE (vxc_lb_h(1:na, 1:nr, 1:nspins), vxc_lb_s(1:na, 1:nr, 1:nspins))
654 ALLOCATE (vxc_gllb_h(1:na, 1:nr, 1:nspins), vxc_gllb_s(1:na, 1:nr, 1:nspins))
655 ALLOCATE (vxc_tmp_h(1:na, 1:nr, 1:nspins), vxc_tmp_s(1:na, 1:nr, 1:nspins))
656 ALLOCATE (vxc_saop_h(1:na, 1:nr, 1:nspins), vxc_saop_s(1:na, 1:nr, 1:nspins))
657 ALLOCATE (drho_h(1:4, 1:na, 1:nr, 1:nspins), drho_s(1:4, 1:na, 1:nr, 1:nspins))
658
659 ! Distribute the atoms of this kind
660 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
661
662 DO ir = 1, nr
663 DO ia = 1, na
664 weight_h(ia, ir) = atomic_grid%wr(ir)*atomic_grid%wa(ia)
665 END DO
666 IF (accint) THEN
667 alpha = dft_control%qs_control%gapw_control%aw(ikind)
668 agr = 1.0_dp - exp(-alpha*atomic_grid%rad2(ir))
669 DO ia = 1, na
670 weight_s(ia, ir) = atomic_grid%wr(ir)*atomic_grid%wa(ia)*agr
671 END DO
672 ELSE
673 DO ia = 1, na
674 weight_s(ia, ir) = atomic_grid%wr(ir)*atomic_grid%wa(ia)
675 END DO
676 END IF
677 END DO
678
679 DO iat = 1, natom !bo(1),bo(2)
680 iatom = atom_list(iat)
681
682 rho_atom => rho_atom_set(iatom)
683 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
684 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
685 rho_rad_s=r_s, drho_rad_h=dr_h, &
686 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
687 rho_rad_s_d=r_s_d)
688 rho_h = 0.0_dp
689 rho_s = 0.0_dp
690 drho_h = 0.0_dp
691 drho_s = 0.0_dp
692 DO ir = 1, nr
693 CALL calc_rho_angular(atomic_grid, harmonics, nspins, .true., &
694 ir, r_h, r_s, rho_h, rho_s, &
695 dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
696 END DO
697 DO ir = 1, nr
698 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau, na, ir)
699 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau, na, ir)
700 END DO
701
702 !-----------------------------!
703 ! 1. Slater exchange for LB94 !
704 !-----------------------------!
705 xc_fun_section_orig => section_vals_get_subs_vals(xc_section_orig, &
706 "XC_FUNCTIONAL")
707 CALL section_vals_create(xc_fun_section_tmp, xc_fun_section_orig%section)
708 CALL section_vals_val_set(xc_fun_section_tmp, "_SECTION_PARAMETERS_", &
709 i_val=xc_funct_no_shortcut)
710 CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
711 l_val=.true.)
712 CALL section_vals_set_subs_vals(xc_section_tmp, "XC_FUNCTIONAL", &
713 xc_fun_section_tmp)
714
715 !---------------------!
716 ! Both: hard and soft !
717 !---------------------!
718 CALL xc_dset_zero_all(deriv_set)
719 CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_h, deriv_set, 1, needs, &
720 weight_h, lsd, na, nr, exc, vxc_tmp_h, vxg, vtau)
721 CALL xc_dset_zero_all(deriv_set)
722 CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_s, deriv_set, 1, needs, &
723 weight_s, lsd, na, nr, exc, vxc_tmp_s, vxg, vtau)
724
725 !-------------------------------------------!
726 ! 2. PZ correlation for LB94 and ec_uniform !
727 !-------------------------------------------!
728 CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
729 l_val=.false.)
730 CALL section_vals_val_set(xc_fun_section_tmp, "PZ81%_SECTION_PARAMETERS_", &
731 l_val=.true.)
732
733 !------!
734 ! Hard !
735 !------!
736 CALL xc_dset_zero_all(deriv_set)
737 CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_h, deriv_set, 1, needs, &
738 weight_h, lsd, na, nr, exc, vxc_lb_h, vxg, vtau)
739 vxc_lb_h = vxc_lb_h + alpha*vxc_tmp_h
740 DO ispin = 1, nspins
741 dummy => vxc_tmp_h(:, :, ispin:ispin)
742 CALL add_lb_pot(dummy, rho_set_h, lsd, ispin)
743 vxc_lb_h(:, :, ispin) = vxc_lb_h(:, :, ispin) - weight_h(:, :)*vxc_tmp_h(:, :, ispin)
744 END DO
745 NULLIFY (dummy)
746
747 vxc_gllb_h = 0.0_dp
748 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
749 cpassert(ASSOCIATED(deriv))
750 CALL xc_derivative_get(deriv, deriv_data=e_uniform)
751 DO ispin = 1, nspins
752 dummy => vxc_gllb_h(:, :, ispin:ispin)
753 CALL calc_2excpbe(dummy, rho_set_h, e_uniform, lsd)
754 vxc_gllb_h(:, :, ispin) = vxc_gllb_h(:, :, ispin)*weight_h(:, :)
755 END DO
756 NULLIFY (deriv, dummy, e_uniform)
757
758 !------!
759 ! Soft !
760 !------!
761 CALL xc_dset_zero_all(deriv_set)
762 CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_s, deriv_set, 1, needs, &
763 weight_s, lsd, na, nr, exc, vxc_lb_s, vxg, vtau)
764
765 vxc_lb_s = vxc_lb_s + alpha*vxc_tmp_s
766 DO ispin = 1, nspins
767 dummy => vxc_tmp_s(:, :, ispin:ispin)
768 CALL add_lb_pot(dummy, rho_set_s, lsd, ispin)
769 vxc_lb_s(:, :, ispin) = vxc_lb_s(:, :, ispin) - weight_s(:, :)*vxc_tmp_s(:, :, ispin)
770 END DO
771 NULLIFY (dummy)
772
773 vxc_gllb_s = 0.0_dp
774 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
775 cpassert(ASSOCIATED(deriv))
776 CALL xc_derivative_get(deriv, deriv_data=e_uniform)
777 DO ispin = 1, nspins
778 dummy => vxc_gllb_s(:, :, ispin:ispin)
779 CALL calc_2excpbe(dummy, rho_set_s, e_uniform, lsd)
780 vxc_gllb_s(:, :, ispin) = vxc_gllb_s(:, :, ispin)*weight_s(:, :)
781 END DO
782 NULLIFY (deriv, dummy, e_uniform)
783
784 !------------------!
785 ! Now the orbitals !
786 !------------------!
787 vxc_tmp_h = 0.0_dp; vxc_tmp_s = 0.0_dp
788
789 DO ispin = 1, nspins
790
791 DO orb = 1, homo(ispin) - 1
792
793 ALLOCATE (coeff_col(nrow(ispin), 1))
794
795 efac = k_rho*sqrt(mo_eigenvalues(ispin)%array(homo(ispin)) - &
796 mo_eigenvalues(ispin)%array(orb))
797 IF (.NOT. lsd) efac = 2.0_dp*efac
798
799 CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
800 CALL cp_fm_get_submatrix(mo_coeff(ispin)%matrix, coeff_col, &
801 1, orb, nrow(ispin), 1)
802 CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
803 1, orb)
804 CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
805 CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
806 matrix_v=single_mo_coeff(ispin), &
807 ncol=ncol(ispin), &
808 alpha=1.0_dp)
809
810 DEALLOCATE (coeff_col)
811
812 ! This calculates the CPC and density on the grids for every atom even though
813 ! we need it only for iatom at the moment. It seems that to circumvent this,
814 ! the routines must be adapted to calculate just iatom
815 ! remap pointer
816 ns = SIZE(orbital_density_matrix)
817 psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
818 CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
819 CALL prepare_gapw_den(qs_env, local_rho_set, .false.)
820
821 rho_atom => local_rho_set%rho_atom_set(iatom)
822 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
823 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
824 rho_h = 0.0_dp
825 rho_s = 0.0_dp
826 drho_h = 0.0_dp
827 drho_s = 0.0_dp
828 DO ir = 1, nr
829 CALL calc_rho_angular(atomic_grid, harmonics, nspins, .false., &
830 ir, r_h, r_s, rho_h, rho_s, &
831 dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
832 END DO
833 DO ir = 1, nr
834 CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
835 CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
836 END DO
837
838 IF (lsd) THEN
839 IF (ispin == 1) THEN
840 vxc_tmp_h(:, :, 1) = vxc_tmp_h(:, :, 1) + efac*orb_rho_set_h%rhoa(:, :, 1)
841 vxc_tmp_s(:, :, 1) = vxc_tmp_s(:, :, 1) + efac*orb_rho_set_s%rhoa(:, :, 1)
842 ELSE
843 vxc_tmp_h(:, :, 2) = vxc_tmp_h(:, :, 2) + efac*orb_rho_set_h%rhob(:, :, 1)
844 vxc_tmp_s(:, :, 2) = vxc_tmp_s(:, :, 2) + efac*orb_rho_set_s%rhob(:, :, 1)
845 END IF
846 ELSE
847 vxc_tmp_h(:, :, 1) = vxc_tmp_h(:, :, 1) + efac*orb_rho_set_h%rho(:, :, 1)
848 vxc_tmp_s(:, :, 1) = vxc_tmp_s(:, :, 1) + efac*orb_rho_set_s%rho(:, :, 1)
849 END IF
850
851 END DO ! orb
852
853 END DO ! ispin
854
855 IF (lsd) THEN
856 DO ir = 1, nr
857 DO ia = 1, na
858 IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
859 vxc_gllb_h(ia, ir, 1) = vxc_gllb_h(ia, ir, 1) + &
860 weight_h(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
861 END IF
862 IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
863 vxc_gllb_h(ia, ir, 2) = vxc_gllb_h(ia, ir, 2) + &
864 weight_h(ia, ir)*vxc_tmp_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
865 END IF
866 IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
867 vxc_gllb_s(ia, ir, 1) = vxc_gllb_s(ia, ir, 1) + &
868 weight_s(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
869 END IF
870 IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
871 vxc_gllb_s(ia, ir, 2) = vxc_gllb_s(ia, ir, 2) + &
872 weight_s(ia, ir)*vxc_tmp_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
873 END IF
874 END DO
875 END DO
876 ELSE
877 DO ir = 1, nr
878 DO ia = 1, na
879 IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
880 vxc_gllb_h(ia, ir, 1) = vxc_gllb_h(ia, ir, 1) + &
881 weight_h(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
882 END IF
883 IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
884 vxc_gllb_s(ia, ir, 1) = vxc_gllb_s(ia, ir, 1) + &
885 weight_s(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
886 END IF
887 END DO
888 END DO
889 END IF
890
891 vxc_saop_h = 0.0_dp; vxc_saop_s = 0.0_dp
892
893 DO ispin = 1, nspins
894
895 DO orb = 1, homo(ispin)
896
897 ALLOCATE (coeff_col(nrow(ispin), 1))
898
899 we_lb = exp(-2.0_dp*(mo_eigenvalues(ispin)%array(homo(ispin)) - &
900 mo_eigenvalues(ispin)%array(orb))**2)
901 we_gllb = 1.0_dp - we_lb
902 IF (.NOT. lsd) THEN
903 we_lb = 2.0_dp*we_lb
904 we_gllb = 2.0_dp*we_gllb
905 END IF
906
907 vxc_tmp_h(:, :, ispin) = we_lb*vxc_lb_h(:, :, ispin) + &
908 we_gllb*vxc_gllb_h(:, :, ispin)
909 vxc_tmp_s(:, :, ispin) = we_lb*vxc_lb_s(:, :, ispin) + &
910 we_gllb*vxc_gllb_s(:, :, ispin)
911
912 CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
913 CALL cp_fm_get_submatrix(mo_coeff(ispin)%matrix, coeff_col, &
914 1, orb, nrow(ispin), 1)
915 CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
916 1, orb)
917 CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
918 CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
919 matrix_v=single_mo_coeff(ispin), &
920 ncol=ncol(ispin), &
921 alpha=1.0_dp)
922
923 DEALLOCATE (coeff_col)
924
925 ! This calculates the CPC and density on the grids for every atom even though
926 ! we need it only for iatom at the moment. It seems that to circumvent this,
927 ! the routines must be adapted to calculate just iatom
928 ! remap pointer
929 ns = SIZE(orbital_density_matrix)
930 psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
931 CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
932 CALL prepare_gapw_den(qs_env, local_rho_set, .false.)
933
934 rho_atom => local_rho_set%rho_atom_set(iatom)
935 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
936 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
937 rho_h = 0.0_dp
938 rho_s = 0.0_dp
939 drho_h = 0.0_dp
940 drho_s = 0.0_dp
941 DO ir = 1, nr
942 CALL calc_rho_angular(atomic_grid, harmonics, nspins, .false., &
943 ir, r_h, r_s, rho_h, rho_s, &
944 dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
945 END DO
946 DO ir = 1, nr
947 CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
948 CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
949 END DO
950
951 IF (lsd) THEN
952 IF (ispin == 1) THEN
953 vxc_saop_h(:, :, 1) = vxc_saop_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rhoa(:, :, 1)
954 vxc_saop_s(:, :, 1) = vxc_saop_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rhoa(:, :, 1)
955 ELSE
956 vxc_saop_h(:, :, 2) = vxc_saop_h(:, :, 2) + vxc_tmp_h(:, :, 2)*orb_rho_set_h%rhob(:, :, 1)
957 vxc_saop_s(:, :, 2) = vxc_saop_s(:, :, 2) + vxc_tmp_s(:, :, 2)*orb_rho_set_s%rhob(:, :, 1)
958 END IF
959 ELSE
960 vxc_saop_h(:, :, 1) = vxc_saop_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rho(:, :, 1)
961 vxc_saop_s(:, :, 1) = vxc_saop_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rho(:, :, 1)
962 END IF
963
964 END DO ! orb
965
966 END DO ! ispin
967
968 IF (lsd) THEN
969 DO ir = 1, nr
970 DO ia = 1, na
971 IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
972 vxc_saop_h(ia, ir, 1) = vxc_saop_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
973 ELSE
974 vxc_saop_h(ia, ir, 1) = 0.0_dp
975 END IF
976 IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
977 vxc_saop_h(ia, ir, 2) = vxc_saop_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
978 ELSE
979 vxc_saop_h(ia, ir, 2) = 0.0_dp
980 END IF
981 IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
982 vxc_saop_s(ia, ir, 1) = vxc_saop_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
983 ELSE
984 vxc_saop_s(ia, ir, 1) = 0.0_dp
985 END IF
986 IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
987 vxc_saop_s(ia, ir, 2) = vxc_saop_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
988 ELSE
989 vxc_saop_s(ia, ir, 2) = 0.0_dp
990 END IF
991 END DO
992 END DO
993 ELSE
994 DO ir = 1, nr
995 DO ia = 1, na
996 IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
997 vxc_saop_h(ia, ir, 1) = vxc_saop_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
998 ELSE
999 vxc_saop_h(ia, ir, 1) = 0.0_dp
1000 END IF
1001 IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
1002 vxc_saop_s(ia, ir, 1) = vxc_saop_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
1003 ELSE
1004 vxc_saop_s(ia, ir, 1) = 0.0_dp
1005 END IF
1006 END DO
1007 END DO
1008 END IF
1009
1010 rho_atom => rho_atom_set(iatom)
1011 CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1012 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis, &
1013 harmonics=harmonics, grid_atom=grid_atom)
1014 SELECT CASE (oe_corr)
1015 CASE (oe_lb)
1016 CALL gavxcgb_nogc(vxc_lb_h, vxc_lb_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1017 CASE (oe_gllb)
1018 CALL gavxcgb_nogc(vxc_gllb_h, vxc_gllb_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1019 CASE (oe_saop)
1020 CALL gavxcgb_nogc(vxc_saop_h, vxc_saop_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1021 CASE default
1022 cpabort("Unknown correction!")
1023 END SELECT
1024
1025 END DO
1026
1027 DEALLOCATE (rho_h, rho_s, weight_h, weight_s)
1028 DEALLOCATE (vxc_lb_h, vxc_lb_s)
1029 DEALLOCATE (vxc_gllb_h, vxc_gllb_s)
1030 DEALLOCATE (vxc_tmp_h, vxc_tmp_s)
1031 DEALLOCATE (vxc_saop_h, vxc_saop_s)
1032 DEALLOCATE (drho_h, drho_s)
1033
1034 CALL xc_dset_release(deriv_set)
1035 CALL xc_rho_set_release(rho_set_h)
1036 CALL xc_rho_set_release(rho_set_s)
1037 CALL xc_rho_set_release(orb_rho_set_h)
1038 CALL xc_rho_set_release(orb_rho_set_s)
1039
1040 END DO
1041
1042 ! remap pointer
1043 ns = SIZE(matrix_ks)
1044 ksmat(1:ns, 1:1) => matrix_ks(1:ns)
1045 ns = SIZE(rho_struct_ao)
1046 psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
1047
1048 CALL update_ks_atom(qs_env, ksmat, psmat, forces=.false.)
1049
1050 !---------!
1051 ! Cleanup !
1052 !---------!
1053 CALL section_vals_release(xc_fun_section_tmp)
1054 CALL section_vals_release(xc_section_tmp)
1055 CALL section_vals_release(xc_section_orig)
1056
1057 CALL local_rho_set_release(local_rho_set)
1058 CALL cp_fm_release(single_mo_coeff)
1059 DEALLOCATE (mo_coeff, mo_eigenvalues)
1060 CALL dbcsr_deallocate_matrix_set(orbital_density_matrix)
1061
1062 END SUBROUTINE gapw_add_atomic_saop_pot
1063
1064! **************************************************************************************************
1065!> \brief ...
1066!> \param pot ...
1067!> \param rho_set ...
1068!> \param lsd ...
1069!> \param spin ...
1070! **************************************************************************************************
1071 SUBROUTINE add_lb_pot(pot, rho_set, lsd, spin)
1072
1073 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pot
1074 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1075 LOGICAL, INTENT(IN) :: lsd
1076 INTEGER, INTENT(IN) :: spin
1077
1078 REAL(kind=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp
1079
1080 INTEGER :: i, j, k
1081 INTEGER, DIMENSION(2, 3) :: bo
1082 REAL(kind=dp) :: n, n_13, x, x2
1083
1084 bo = rho_set%local_bounds
1085
1086 DO k = bo(1, 3), bo(2, 3)
1087 DO j = bo(1, 2), bo(2, 2)
1088 DO i = bo(1, 1), bo(2, 1)
1089 IF (.NOT. lsd) THEN
1090 IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1091 n = rho_set%rho(i, j, k)/2.0_dp
1092 n_13 = n**ob3
1093 x = (rho_set%norm_drho(i, j, k)/2.0_dp)/(n*n_13)
1094 x2 = x*x
1095 pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(x + sqrt(x2 + 1.0_dp)))
1096 END IF
1097 ELSE
1098 IF (spin == 1) THEN
1099 IF (rho_set%rhoa(i, j, k) > rho_set%rho_cutoff) THEN
1100 n_13 = rho_set%rhoa_1_3(i, j, k)
1101 x = rho_set%norm_drhoa(i, j, k)/(rho_set%rhoa(i, j, k)*n_13)
1102 x2 = x*x
1103 pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(sqrt(x2 + 1.0_dp) + x))
1104 END IF
1105 ELSE IF (spin == 2) THEN
1106 IF (rho_set%rhob(i, j, k) > rho_set%rho_cutoff) THEN
1107 n_13 = rho_set%rhob_1_3(i, j, k)
1108 x = rho_set%norm_drhob(i, j, k)/(rho_set%rhob(i, j, k)*n_13)
1109 x2 = x*x
1110 pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(sqrt(x2 + 1.0_dp) + x))
1111 END IF
1112 END IF
1113 END IF
1114 END DO
1115 END DO
1116 END DO
1117
1118 END SUBROUTINE add_lb_pot
1119
1120! **************************************************************************************************
1121!> \brief ...
1122!> \param pot ...
1123!> \param rho_set ...
1124!> \param e_uniform ...
1125!> \param lsd ...
1126! **************************************************************************************************
1127 SUBROUTINE calc_2excpbe(pot, rho_set, e_uniform, lsd)
1128
1129 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pot
1130 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1131 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: e_uniform
1132 LOGICAL, INTENT(IN) :: lsd
1133
1134 INTEGER :: i, j, k
1135 INTEGER, DIMENSION(2, 3) :: bo
1136 REAL(kind=dp) :: e_unif, rho
1137
1138 bo = rho_set%local_bounds
1139
1140 DO k = bo(1, 3), bo(2, 3)
1141 DO j = bo(1, 2), bo(2, 2)
1142 DO i = bo(1, 1), bo(2, 1)
1143 IF (.NOT. lsd) THEN
1144 IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1145 e_unif = e_uniform(i, j, k)/rho_set%rho(i, j, k)
1146 ELSE
1147 e_unif = 0.0_dp
1148 END IF
1149 pot(i, j, k) = &
1150 2.0_dp* &
1151 calc_ecpbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1152 e_unif, rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1153 2.0_dp* &
1154 calc_expbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1155 rho_set%rho_cutoff, rho_set%drho_cutoff)
1156 ELSE
1157 rho = rho_set%rhoa(i, j, k) + rho_set%rhob(i, j, k)
1158 IF (rho > rho_set%rho_cutoff) THEN
1159 e_unif = e_uniform(i, j, k)/rho
1160 ELSE
1161 e_unif = 0.0_dp
1162 END IF
1163 pot(i, j, k) = &
1164 2.0_dp* &
1165 calc_ecpbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1166 e_unif, &
1167 rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1168 2.0_dp* &
1169 calc_expbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1170 rho_set%rho_cutoff, rho_set%drho_cutoff)
1171 END IF
1172 END DO
1173 END DO
1174 END DO
1175
1176 END SUBROUTINE calc_2excpbe
1177
1178! **************************************************************************************************
1179!> \brief ...
1180!> \param ra ...
1181!> \param rb ...
1182!> \param ngr ...
1183!> \param ec_unif ...
1184!> \param rc ...
1185!> \param ngrc ...
1186!> \return ...
1187! **************************************************************************************************
1188 FUNCTION calc_ecpbe_u(ra, rb, ngr, ec_unif, rc, ngrc) RESULT(res)
1189
1190 REAL(kind=dp), INTENT(in) :: ra, rb, ngr, ec_unif, rc, ngrc
1191 REAL(kind=dp) :: res
1192
1193 REAL(kind=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp, tb3 = 2.0_dp/3.0_dp
1194
1195 REAL(kind=dp) :: a, at2, h, kf, kl, ks, phi, phi3, r, t2, &
1196 zeta
1197
1198 r = ra + rb
1199 h = 0.0_dp
1200 IF (r > rc .AND. ngr > ngrc) THEN
1201 zeta = (ra - rb)/r
1202 IF (zeta > 1.0_dp) zeta = 1.0_dp ! machine precision problem
1203 IF (zeta < -1.0_dp) zeta = -1.0_dp ! machine precision problem
1204 phi = ((1.0_dp + zeta)**tb3 + (1.0_dp - zeta)**tb3)/2.0_dp
1205 phi3 = phi*phi*phi
1206 kf = (3.0_dp*r*pi*pi)**ob3
1207 ks = sqrt(4.0_dp*kf/pi)
1208 t2 = (ngr/(2.0_dp*phi*ks*r))**2
1209 a = beta_ec/gamma_saop/(exp(-ec_unif/(gamma_saop*phi3)) - 1.0_dp)
1210 at2 = a*t2
1211 kl = (1.0_dp + at2)/(1.0_dp + at2 + at2*at2)
1212 h = gamma_saop*log(1.0_dp + beta_ec/gamma_saop*t2*kl)
1213 END IF
1214 res = ec_unif + h
1215
1216 END FUNCTION calc_ecpbe_u
1217
1218! **************************************************************************************************
1219!> \brief ...
1220!> \param r ...
1221!> \param ngr ...
1222!> \param ec_unif ...
1223!> \param rc ...
1224!> \param ngrc ...
1225!> \return ...
1226! **************************************************************************************************
1227 FUNCTION calc_ecpbe_r(r, ngr, ec_unif, rc, ngrc) RESULT(res)
1228
1229 REAL(kind=dp), INTENT(in) :: r, ngr, ec_unif, rc, ngrc
1230 REAL(kind=dp) :: res
1231
1232 REAL(kind=dp) :: a, at2, h, kf, kl, ks, t2
1233
1234 h = 0.0_dp
1235 IF (r > rc .AND. ngr > ngrc) THEN
1236 kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1237 ks = sqrt(4.0_dp*kf/pi)
1238 t2 = (ngr/(2.0_dp*ks*r))**2
1239 a = beta_ec/gamma_saop/(exp(-ec_unif/gamma_saop) - 1.0_dp)
1240 at2 = a*t2
1241 kl = (1.0_dp + at2)/(1.0_dp + at2 + at2*at2)
1242 h = gamma_saop*log(1.0_dp + beta_ec/gamma_saop*t2*kl)
1243 END IF
1244 res = ec_unif + h
1245
1246 END FUNCTION calc_ecpbe_r
1247
1248! **************************************************************************************************
1249!> \brief ...
1250!> \param ra ...
1251!> \param rb ...
1252!> \param ngr ...
1253!> \param rc ...
1254!> \param ngrc ...
1255!> \return ...
1256! **************************************************************************************************
1257 FUNCTION calc_expbe_u(ra, rb, ngr, rc, ngrc) RESULT(res)
1258
1259 REAL(kind=dp), INTENT(in) :: ra, rb, ngr, rc, ngrc
1260 REAL(kind=dp) :: res
1261
1262 REAL(kind=dp) :: r
1263
1264 r = ra + rb
1265 res = calc_expbe_r(r, ngr, rc, ngrc)
1266
1267 END FUNCTION calc_expbe_u
1268
1269! **************************************************************************************************
1270!> \brief ...
1271!> \param r ...
1272!> \param ngr ...
1273!> \param rc ...
1274!> \param ngrc ...
1275!> \return ...
1276! **************************************************************************************************
1277 FUNCTION calc_expbe_r(r, ngr, rc, ngrc) RESULT(res)
1278
1279 REAL(kind=dp), INTENT(in) :: r, ngr, rc, ngrc
1280 REAL(kind=dp) :: res
1281
1282 REAL(kind=dp) :: ex_unif, fx, kf, s
1283
1284 IF (r > rc) THEN
1285 kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1286 ex_unif = -3.0_dp*kf/(4.0_dp*pi)
1287 fx = 1.0_dp
1288 IF (ngr > ngrc) THEN
1289 s = ngr/(2.0_dp*kf*r)
1290 fx = fx + kappa - kappa/(1.0_dp + mu*s*s/kappa)
1291 END IF
1292 res = ex_unif*fx
1293 ELSE
1294 res = 0.0_dp
1295 END IF
1296
1297 END FUNCTION calc_expbe_r
1298
1299END MODULE xc_pot_saop
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:88
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_zero
Definition pw_methods.F:83
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_deallocate_matrix
subroutine, public dbcsr_deallocate_matrix(matrix)
...
Definition cp_dbcsr_api.F:235
cp_dbcsr_api::dbcsr_copy
subroutine, public dbcsr_copy(matrix_b, matrix_a, name, keep_sparsity, keep_imaginary)
...
Definition cp_dbcsr_api.F:371
cp_dbcsr_api::dbcsr_set
subroutine, public dbcsr_set(matrix, alpha)
...
Definition cp_dbcsr_api.F:1204
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_dbcsr_operations::cp_dbcsr_plus_fm_fm_t
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,...
Definition cp_dbcsr_operations.F:693
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_get_info
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
Definition cp_fm_types.F:1038
cp_fm_types::cp_fm_set_submatrix
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition cp_fm_types.F:750
cp_fm_types::cp_fm_set_all
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
Definition cp_fm_types.F:534
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, nrow, ncol, set_zero)
creates a new full matrix with the given structure
Definition cp_fm_types.F:165
cp_fm_types::cp_fm_get_submatrix
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition cp_fm_types.F:925
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::oe_saop
integer, parameter, public oe_saop
Definition input_constants.F:608
input_constants::xc_funct_no_shortcut
integer, parameter, public xc_funct_no_shortcut
Definition input_constants.F:572
input_constants::do_method_gapw
integer, parameter, public do_method_gapw
Definition input_constants.F:260
input_constants::oe_lb
integer, parameter, public oe_lb
Definition input_constants.F:608
input_constants::oe_gllb
integer, parameter, public oe_gllb
Definition input_constants.F:608
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_create
recursive subroutine, public section_vals_create(section_vals, section)
creates a object where to store the values of a section
Definition input_section_types.F:612
input_section_types::section_vals_val_set
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
Definition input_section_types.F:1233
input_section_types::section_vals_retain
subroutine, public section_vals_retain(section_vals)
retains the given section values (see doc/ReferenceCounting.html)
Definition input_section_types.F:642
input_section_types::section_vals_get_subs_vals
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
Definition input_section_types.F:735
input_section_types::section_vals_set_subs_vals
subroutine, public section_vals_set_subs_vals(section_vals, subsection_name, new_section_vals, i_rep_section)
replaces of the requested subsection with the one given
Definition input_section_types.F:1723
input_section_types::section_vals_duplicate
subroutine, public section_vals_duplicate(section_vals_in, section_vals_out, i_rep_start, i_rep_end)
creates a deep copy from section_vals_in to section_vals_out
Definition input_section_types.F:1777
input_section_types::section_vals_val_get
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
Definition input_section_types.F:1052
input_section_types::section_vals_release
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Definition input_section_types.F:657
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
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
Definition pw_env_types.F:113
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::calculate_rho_elec
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
Definition qs_collocate_density.F:1720
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
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, mimic, 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, xcint_weights, 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.
Definition qs_environment_types.F:533
qs_gapw_densities
Definition qs_gapw_densities.F:9
qs_gapw_densities::prepare_gapw_den
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external, pw_env_sub)
...
Definition qs_gapw_densities.F:57
qs_grid_atom
Definition qs_grid_atom.F:8
qs_harmonics_atom
Definition qs_harmonics_atom.F:8
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
qs_ks_atom
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
Definition qs_ks_atom.F:12
qs_ks_atom::update_ks_atom
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
Definition qs_ks_atom.F:110
qs_ks_types
Definition qs_ks_types.F:15
qs_local_rho_types
Definition qs_local_rho_types.F:8
qs_local_rho_types::local_rho_set_create
subroutine, public local_rho_set_create(local_rho_set)
...
Definition qs_local_rho_types.F:203
qs_local_rho_types::local_rho_set_release
subroutine, public local_rho_set_release(local_rho_set)
...
Definition qs_local_rho_types.F:224
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::get_mo_set
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.
Definition qs_mo_types.F:401
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_oce_types
Definition qs_oce_types.F:8
qs_rho_atom_methods
Definition qs_rho_atom_methods.F:7
qs_rho_atom_methods::allocate_rho_atom_internals
subroutine, public allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
Definition qs_rho_atom_methods.F:960
qs_rho_atom_methods::calculate_rho_atom_coeff
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
Definition qs_rho_atom_methods.F:454
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_atom_types::get_rho_atom
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_vlocal_gb_h, ga_vlocal_gb_s, int_scr_h, int_scr_s)
...
Definition qs_rho_atom_types.F:195
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
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...
Definition qs_rho_types.F:229
qs_vxc_atom
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition qs_vxc_atom.F:12
qs_vxc_atom::calc_rho_angular
subroutine, public calc_rho_angular(grid_atom, harmonics, nspins, grad_func, ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
...
Definition qs_vxc_atom.F:2080
qs_vxc_atom::gavxcgb_nogc
subroutine, public gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition qs_vxc_atom.F:2445
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
virial_types
Definition virial_types.F:13
xc_atom
Definition xc_atom.F:9
xc_atom::vxc_of_r_new
subroutine, public vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, lsd, na, nr, exc, vxc, vxg, vtau, energy_only, adiabatic_rescale_factor)
...
Definition xc_atom.F:64
xc_atom::xc_rho_set_atom_update
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition xc_atom.F:507
xc_atom::fill_rho_set
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition xc_atom.F:667
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_zero_all
subroutine, public xc_dset_zero_all(deriv_set)
...
Definition xc_derivative_set_types.F:166
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:65
xc_derivative_set_types::xc_dset_release
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
Definition xc_derivative_set_types.F:143
xc_derivative_set_types::xc_dset_create
subroutine, public xc_dset_create(derivative_set, pw_pool, local_bounds)
creates a derivative set object
Definition xc_derivative_set_types.F:111
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_derivatives
Definition xc_derivatives.F:9
xc_derivatives::xc_functionals_eval
subroutine, public xc_functionals_eval(functionals, lsd, rho_set, deriv_set, deriv_order)
...
Definition xc_derivatives.F:735
xc_pot_saop
Calculate the saop potential.
Definition xc_pot_saop.F:11
xc_pot_saop::add_saop_pot
subroutine, public add_saop_pot(ks_matrix, qs_env, oe_corr)
...
Definition xc_pot_saop.F:126
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_cflags_types::xc_rho_cflags_setall
elemental subroutine, public xc_rho_cflags_setall(cflags, value)
sets all the flags to the given value
Definition xc_rho_cflags_types.F:63
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_create
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
Definition xc_rho_set_types.F:111
xc_rho_set_types::xc_rho_set_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition xc_rho_set_types.F:129
xc_rho_set_types::xc_rho_set_update
subroutine, public xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, xc_deriv_method_id, xc_rho_smooth_id, pw_pool, spinflip)
updates the given rho set with the density given by rho_r (and rho_g). The rho set will contain the c...
Definition xc_rho_set_types.F:691
xc_xbecke88
calculates the Becke 88 exchange functional
Definition xc_xbecke88.F:14
xc_xbecke88::xb88_lsd_info
subroutine, public xb88_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_xbecke88.F:90
xc_xbecke88::xb88_lda_info
subroutine, public xb88_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_xbecke88.F:59
xc
Exchange and Correlation functional calculations.
Definition xc.F:17
xc::xc_vxc_pw_create
subroutine, public xc_vxc_pw_create(vxc_rho, vxc_tau, exc, rho_r, rho_g, tau, xc_section, weights, pw_pool, compute_virial, virial_xc, exc_r)
Exchange and Correlation functional calculations.
Definition xc.F:474
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_array_utils::cp_1d_r_p_type
represent a pointer to a 1d array
Definition cp_array_utils.F:99
cp_control_types::dft_control_type
Definition cp_control_types.F:751
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
cp_fm_types::cp_fm_p_type
just to build arrays of pointers to matrices
Definition cp_fm_types.F:128
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:114
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:746
orbital
Orbital angular momentum.
Definition grid_common.h:133
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44
qs_harmonics_atom::harmonics_atom_type
Definition qs_harmonics_atom.F:32
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_local_rho_types::local_rho_type
Definition qs_local_rho_types.F:45
qs_mo_types::mo_set_type
Definition qs_mo_types.F:47
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_oce_types::oce_matrix_type
Definition qs_oce_types.F:37
qs_rho_atom_types::rho_atom_coeff
Definition qs_rho_atom_types.F:19
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
virial_types::virial_type
Definition virial_types.F:26
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78