(git:6344619)
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_utils, ONLY: calc_rho_angular,&
79 calc_weight_function,&
80 gavxcgb_nogc
81 USE util, ONLY: get_limit
82 USE virial_types, ONLY: virial_type
83 USE xc, ONLY: xc_vxc_pw_create
84 USE xc_atom, ONLY: fill_rho_set,&
85 vxc_of_r_new,&
86 xc_rho_set_atom_update
87 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
88 xc_dset_create,&
89 xc_dset_get_derivative,&
90 xc_dset_release,&
91 xc_dset_zero_all
92 USE xc_derivative_types, ONLY: xc_derivative_get,&
93 xc_derivative_type
94 USE xc_derivatives, ONLY: xc_functionals_eval
95 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_setall,&
96 xc_rho_cflags_type
97 USE xc_rho_set_types, ONLY: xc_rho_set_create,&
98 xc_rho_set_release,&
99 xc_rho_set_type,&
100 xc_rho_set_update
101 USE xc_xbecke88, ONLY: xb88_lda_info,&
102 xb88_lsd_info
103#include "./base/base_uses.f90"
104
105 IMPLICIT NONE
106
107 PRIVATE
108
109 PUBLIC :: add_saop_pot
110
111 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pot_saop'
112
113 ! should be eliminated
114 REAL(KIND=dp), PARAMETER :: alpha = 1.19_dp, beta = 0.01_dp, k_rho = 0.42_dp
115 REAL(KIND=dp), PARAMETER :: kappa = 0.804_dp, mu = 0.21951_dp, &
116 beta_ec = 0.066725_dp, gamma_saop = 0.031091_dp
117
118CONTAINS
119
120! **************************************************************************************************
121!> \brief ...
122!> \param ks_matrix ...
123!> \param qs_env ...
124!> \param oe_corr ...
125! **************************************************************************************************
126 SUBROUTINE add_saop_pot(ks_matrix, qs_env, oe_corr)
127
128 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
129 TYPE(qs_environment_type), POINTER :: qs_env
130 INTEGER, INTENT(IN) :: oe_corr
131
132 INTEGER :: dft_method_id, homo, i, ispin, j, k, &
133 nspins, orb, xc_deriv_method_id, &
134 xc_rho_smooth_id
135 INTEGER, DIMENSION(2) :: ncol, nrow
136 INTEGER, DIMENSION(2, 3) :: bo
137 LOGICAL :: compute_virial, gapw, lsd
138 REAL(kind=dp) :: density_cut, efac, gradient_cut, &
139 tau_cut, we_gllb, we_lb, xc_energy
140 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_col
141 REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
142 REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
143 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: e_uniform
144 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: single_mo_coeff
145 TYPE(cp_fm_type), POINTER :: mo_coeff
146 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: orbital_density_matrix, rho_struct_ao
147 TYPE(mo_set_type), DIMENSION(:), POINTER :: molecular_orbitals
148 TYPE(pw_c1d_gs_type) :: orbital_g
149 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
150 TYPE(pw_env_type), POINTER :: pw_env
151 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
152 TYPE(pw_r3d_rs_type) :: orbital
153 TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: vxc_gllb, vxc_saop
154 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_struct_r, tau, vxc_lb, &
155 vxc_tau, vxc_tmp
156 TYPE(pw_r3d_rs_type), POINTER :: weights
157 TYPE(qs_ks_env_type), POINTER :: ks_env
158 TYPE(qs_rho_type), POINTER :: rho_struct
159 TYPE(section_vals_type), POINTER :: input, xc_fun_section_orig, &
160 xc_fun_section_tmp, xc_section_orig, &
161 xc_section_tmp
162 TYPE(virial_type), POINTER :: virial
163 TYPE(xc_derivative_set_type) :: deriv_set
164 TYPE(xc_derivative_type), POINTER :: deriv
165 TYPE(xc_rho_cflags_type) :: needs
166 TYPE(xc_rho_set_type) :: rho_set
167
168 NULLIFY (ks_env, pw_env, auxbas_pw_pool, input)
169 NULLIFY (rho_g, rho_r, tau, rho_struct, e_uniform)
170 NULLIFY (vxc_lb, vxc_tmp, vxc_tau)
171 NULLIFY (mo_eigenvalues, deriv, rho_struct_r, rho_struct_ao)
172 NULLIFY (orbital_density_matrix, xc_section_tmp, xc_fun_section_tmp)
173
174 CALL get_qs_env(qs_env, &
175 ks_env=ks_env, &
176 rho=rho_struct, &
177 xcint_weights=weights, &
178 pw_env=pw_env, &
179 input=input, &
180 virial=virial, &
181 mos=molecular_orbitals)
182 compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
183 CALL section_vals_val_get(input, "DFT%QS%METHOD", i_val=dft_method_id)
184 gapw = (dft_method_id == do_method_gapw)
185
186 xc_section_orig => section_vals_get_subs_vals(input, "DFT%XC")
187 CALL section_vals_retain(xc_section_orig)
188 CALL section_vals_duplicate(xc_section_orig, xc_section_tmp)
189
190 CALL section_vals_val_get(xc_section_orig, "DENSITY_CUTOFF", &
191 r_val=density_cut)
192 CALL section_vals_val_get(xc_section_orig, "GRADIENT_CUTOFF", &
193 r_val=gradient_cut)
194 CALL section_vals_val_get(xc_section_orig, "TAU_CUTOFF", &
195 r_val=tau_cut)
196
197 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
198
199 CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
200 IF (lsd) THEN
201 nspins = 2
202 ELSE
203 nspins = 1
204 END IF
205
206 ALLOCATE (single_mo_coeff(nspins))
207 CALL dbcsr_allocate_matrix_set(orbital_density_matrix, nspins)
208 CALL qs_rho_get(rho_struct, rho_r=rho_struct_r, rho_ao=rho_struct_ao)
209 rho_r => rho_struct_r
210 DO ispin = 1, nspins
211 ALLOCATE (orbital_density_matrix(ispin)%matrix)
212 CALL dbcsr_copy(orbital_density_matrix(ispin)%matrix, &
213 rho_struct_ao(ispin)%matrix, "orbital density")
214 END DO
215 bo = rho_r(1)%pw_grid%bounds_local
216
217 !---------------------------!
218 ! create the density needed !
219 !---------------------------!
220 CALL xc_rho_set_create(rho_set, bo, &
221 density_cut, &
222 gradient_cut, &
223 tau_cut)
224 CALL xc_rho_cflags_setall(needs, .false.)
225 IF (lsd) THEN
226 CALL xb88_lsd_info(needs=needs)
227 needs%norm_drho = .true.
228 ELSE
229 CALL xb88_lda_info(needs=needs)
230 END IF
231 CALL section_vals_val_get(xc_section_orig, "XC_GRID%XC_DERIV", &
232 i_val=xc_deriv_method_id)
233 CALL section_vals_val_get(xc_section_orig, "XC_GRID%XC_SMOOTH_RHO", &
234 i_val=xc_rho_smooth_id)
235 CALL xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, &
236 xc_deriv_method_id, &
237 xc_rho_smooth_id, &
238 auxbas_pw_pool)
239
240 !----------------------------------------!
241 ! Construct the LB94 potential in vxc_LB !
242 !----------------------------------------!
243 xc_fun_section_orig => section_vals_get_subs_vals(xc_section_orig, &
244 "XC_FUNCTIONAL")
245 CALL section_vals_create(xc_fun_section_tmp, xc_fun_section_orig%section)
246 CALL section_vals_val_set(xc_fun_section_tmp, "_SECTION_PARAMETERS_", &
247 i_val=xc_funct_no_shortcut)
248 CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
249 l_val=.true.)
250 CALL section_vals_set_subs_vals(xc_section_tmp, "XC_FUNCTIONAL", &
251 xc_fun_section_tmp)
252
253 cpassert(.NOT. compute_virial)
254 CALL xc_vxc_pw_create(vxc_tmp, vxc_tau, xc_energy, rho_r, rho_g, tau, &
255 xc_section_tmp, weights, auxbas_pw_pool, &
256 compute_virial=.false., virial_xc=virial_xc_tmp)
257
258 CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
259 l_val=.false.)
260 CALL section_vals_val_set(xc_fun_section_tmp, "PZ81%_SECTION_PARAMETERS_", &
261 l_val=.true.)
262
263 cpassert(.NOT. compute_virial)
264 CALL xc_vxc_pw_create(vxc_lb, vxc_tau, xc_energy, rho_r, rho_g, tau, &
265 xc_section_tmp, weights, auxbas_pw_pool, &
266 compute_virial=.false., virial_xc=virial_xc_tmp)
267
268 DO ispin = 1, nspins
269 CALL pw_axpy(vxc_tmp(ispin), vxc_lb(ispin), alpha)
270 END DO
271
272 DO ispin = 1, nspins
273 CALL add_lb_pot(vxc_tmp(ispin)%array, rho_set, lsd, ispin)
274 CALL pw_axpy(vxc_tmp(ispin), vxc_lb(ispin), -1.0_dp)
275 END DO
276
277 !-----------------------------------------------------------------------------------!
278 ! Construct 2 times PBE one particle density from the PZ correlation energy density !
279 !-----------------------------------------------------------------------------------!
280 CALL xc_dset_create(deriv_set, local_bounds=bo)
281 CALL xc_functionals_eval(xc_fun_section_tmp, &
282 lsd=lsd, &
283 rho_set=rho_set, &
284 deriv_set=deriv_set, &
285 deriv_order=0)
286
287 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
288 CALL xc_derivative_get(deriv, deriv_data=e_uniform)
289
290 ALLOCATE (vxc_gllb(nspins))
291 DO ispin = 1, nspins
292 CALL auxbas_pw_pool%create_pw(vxc_gllb(ispin))
293 END DO
294
295 DO ispin = 1, nspins
296 CALL calc_2excpbe(vxc_gllb(ispin)%array, rho_set, e_uniform, lsd)
297 END DO
298
299 CALL xc_dset_release(deriv_set)
300
301 CALL auxbas_pw_pool%create_pw(orbital)
302 CALL auxbas_pw_pool%create_pw(orbital_g)
303
304 DO ispin = 1, nspins
305
306 CALL get_mo_set(molecular_orbitals(ispin), &
307 mo_coeff=mo_coeff, &
308 eigenvalues=mo_eigenvalues, &
309 homo=homo)
310 CALL cp_fm_create(single_mo_coeff(ispin), &
311 mo_coeff%matrix_struct, &
312 "orbital density matrix")
313
314 CALL cp_fm_get_info(single_mo_coeff(ispin), &
315 nrow_global=nrow(ispin), ncol_global=ncol(ispin))
316 ALLOCATE (coeff_col(nrow(ispin), 1))
317
318 CALL pw_zero(vxc_tmp(ispin))
319
320 DO orb = 1, homo - 1
321
322 efac = k_rho*sqrt(mo_eigenvalues(homo) - mo_eigenvalues(orb))
323 IF (.NOT. lsd) efac = 2.0_dp*efac
324
325 CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
326 CALL cp_fm_get_submatrix(mo_coeff, coeff_col, &
327 1, orb, nrow(ispin), 1)
328 CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
329 1, orb)
330 CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
331 CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
332 matrix_v=single_mo_coeff(ispin), &
333 ncol=ncol(ispin), &
334 alpha=1.0_dp)
335 CALL pw_zero(orbital)
336 CALL pw_zero(orbital_g)
337 CALL calculate_rho_elec(matrix_p=orbital_density_matrix(ispin)%matrix, &
338 rho=orbital, rho_gspace=orbital_g, &
339 ks_env=ks_env)
340
341 CALL pw_axpy(orbital, vxc_tmp(ispin), efac)
342
343 END DO
344 DEALLOCATE (coeff_col)
345
346 DO k = bo(1, 3), bo(2, 3)
347 DO j = bo(1, 2), bo(2, 2)
348 DO i = bo(1, 1), bo(2, 1)
349 IF (rho_r(ispin)%array(i, j, k) > density_cut) THEN
350 vxc_tmp(ispin)%array(i, j, k) = vxc_tmp(ispin)%array(i, j, k)/ &
351 rho_r(ispin)%array(i, j, k)
352 ELSE
353 vxc_tmp(ispin)%array(i, j, k) = 0.0_dp
354 END IF
355 END DO
356 END DO
357 END DO
358
359 CALL pw_axpy(vxc_tmp(ispin), vxc_gllb(ispin), 1.0_dp)
360
361 END DO
362
363 !---------------!
364 ! Assemble SAOP !
365 !---------------!
366 ALLOCATE (vxc_saop(nspins))
367
368 DO ispin = 1, nspins
369
370 CALL get_mo_set(molecular_orbitals(ispin), &
371 mo_coeff=mo_coeff, &
372 eigenvalues=mo_eigenvalues, &
373 homo=homo)
374 CALL auxbas_pw_pool%create_pw(vxc_saop(ispin))
375 CALL pw_zero(vxc_saop(ispin))
376
377 ALLOCATE (coeff_col(nrow(ispin), 1))
378
379 DO orb = 1, homo
380
381 we_lb = exp(-2.0_dp*(mo_eigenvalues(homo) - mo_eigenvalues(orb))**2)
382 we_gllb = 1.0_dp - we_lb
383 IF (.NOT. lsd) THEN
384 we_lb = 2.0_dp*we_lb
385 we_gllb = 2.0_dp*we_gllb
386 END IF
387
388 vxc_tmp(ispin)%array = we_lb*vxc_lb(ispin)%array + &
389 we_gllb*vxc_gllb(ispin)%array
390
391 CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
392 CALL cp_fm_get_submatrix(mo_coeff, coeff_col, &
393 1, orb, nrow(ispin), 1)
394 CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
395 1, orb)
396 CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
397 CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
398 matrix_v=single_mo_coeff(ispin), &
399 ncol=ncol(ispin), &
400 alpha=1.0_dp)
401 CALL pw_zero(orbital)
402 CALL pw_zero(orbital_g)
403 CALL calculate_rho_elec(matrix_p=orbital_density_matrix(ispin)%matrix, &
404 rho=orbital, rho_gspace=orbital_g, &
405 ks_env=ks_env)
406
407 vxc_saop(ispin)%array = vxc_saop(ispin)%array + &
408 orbital%array*vxc_tmp(ispin)%array
409
410 END DO
411
412 CALL dbcsr_deallocate_matrix(orbital_density_matrix(ispin)%matrix)
413
414 DEALLOCATE (coeff_col)
415
416 DO k = bo(1, 3), bo(2, 3)
417 DO j = bo(1, 2), bo(2, 2)
418 DO i = bo(1, 1), bo(2, 1)
419 IF (rho_r(ispin)%array(i, j, k) > density_cut) THEN
420 vxc_saop(ispin)%array(i, j, k) = vxc_saop(ispin)%array(i, j, k)/ &
421 rho_r(ispin)%array(i, j, k)
422 ELSE
423 vxc_saop(ispin)%array(i, j, k) = 0.0_dp
424 END IF
425 END DO
426 END DO
427 END DO
428
429 END DO
430
431 CALL cp_fm_release(single_mo_coeff)
432
433 CALL xc_rho_set_release(rho_set, auxbas_pw_pool)
434 CALL auxbas_pw_pool%give_back_pw(orbital)
435 CALL auxbas_pw_pool%give_back_pw(orbital_g)
436
437 !--------------------!
438 ! Do the integration !
439 !--------------------!
440 DO ispin = 1, nspins
441
442 IF (oe_corr == oe_lb) THEN
443 CALL pw_copy(vxc_lb(ispin), vxc_saop(ispin))
444 ELSE IF (oe_corr == oe_gllb) THEN
445 CALL pw_copy(vxc_gllb(ispin), vxc_saop(ispin))
446 END IF
447 CALL pw_scale(vxc_saop(ispin), vxc_saop(ispin)%pw_grid%dvol)
448
449 CALL integrate_v_rspace(v_rspace=vxc_saop(ispin), pmat=rho_struct_ao(ispin), &
450 hmat=ks_matrix(ispin), qs_env=qs_env, &
451 calculate_forces=.false., &
452 gapw=gapw)
453
454 END DO
455
456 DO ispin = 1, nspins
457 CALL auxbas_pw_pool%give_back_pw(vxc_saop(ispin))
458 CALL auxbas_pw_pool%give_back_pw(vxc_gllb(ispin))
459 CALL vxc_lb(ispin)%release()
460 CALL vxc_tmp(ispin)%release()
461 END DO
462 DEALLOCATE (vxc_gllb, vxc_lb, vxc_tmp, orbital_density_matrix)
463
464 DEALLOCATE (vxc_saop)
465
466 CALL section_vals_release(xc_fun_section_tmp)
467 CALL section_vals_release(xc_section_tmp)
468 CALL section_vals_release(xc_section_orig)
469
470 !-----------------------!
471 ! Call the GAPW routine !
472 !-----------------------!
473 IF (gapw) THEN
474 CALL gapw_add_atomic_saop_pot(qs_env, oe_corr)
475 END IF
476
477 END SUBROUTINE add_saop_pot
478
479! **************************************************************************************************
480!> \brief ...
481!> \param qs_env ...
482!> \param oe_corr ...
483! **************************************************************************************************
484 SUBROUTINE gapw_add_atomic_saop_pot(qs_env, oe_corr)
485
486 TYPE(qs_environment_type), POINTER :: qs_env
487 INTEGER, INTENT(IN) :: oe_corr
488
489 INTEGER :: ia, iat, iatom, ikind, ir, ispin, na, &
490 natom, nr, ns, nspins, on, orb
491 INTEGER, DIMENSION(2) :: bo, homo, ncol, nrow
492 INTEGER, DIMENSION(2, 3) :: bounds
493 INTEGER, DIMENSION(:), POINTER :: atom_list
494 LOGICAL :: accint, lsd, paw_atom
495 REAL(dp), DIMENSION(:, :, :), POINTER :: tau
496 REAL(kind=dp) :: agr, alpha, density_cut, efac, exc, &
497 gradient_cut, tau_cut, we_gllb, we_lb
498 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: fw
499 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_col, weight_h, weight_s
500 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_uniform, rho_h, rho_s, vtau, &
501 vxc_gllb_h, vxc_gllb_s, vxc_lb_h, vxc_lb_s, vxc_saop_h, vxc_saop_s, vxc_tmp_h, vxc_tmp_s
502 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg
503 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
504 TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: mo_eigenvalues
505 TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff
506 TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: single_mo_coeff
507 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, orbital_density_matrix, &
508 rho_struct_ao
509 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat
510 TYPE(dft_control_type), POINTER :: dft_control
511 TYPE(grid_atom_type), POINTER :: atomic_grid, grid_atom
512 TYPE(gto_basis_set_type), POINTER :: orb_basis
513 TYPE(harmonics_atom_type), POINTER :: harmonics
514 TYPE(local_rho_type), POINTER :: local_rho_set
515 TYPE(mo_set_type), DIMENSION(:), POINTER :: molecular_orbitals
516 TYPE(mp_para_env_type), POINTER :: para_env
517 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
518 POINTER :: sab
519 TYPE(oce_matrix_type), POINTER :: oce
520 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
521 TYPE(qs_rho_type), POINTER :: rho_structure
522 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
523 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
524 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
525 TYPE(rho_atom_type), POINTER :: rho_atom
526 TYPE(section_vals_type), POINTER :: input, xc_fun_section_orig, &
527 xc_fun_section_tmp, xc_section_orig, &
528 xc_section_tmp
529 TYPE(xc_derivative_set_type) :: deriv_set
530 TYPE(xc_derivative_type), POINTER :: deriv
531 TYPE(xc_rho_cflags_type) :: needs, needs_orbs
532 TYPE(xc_rho_set_type) :: orb_rho_set_h, orb_rho_set_s, rho_set_h, &
533 rho_set_s
534
535 NULLIFY (rho_h, rho_s, vxc_lb_h, vxc_lb_s, vxc_gllb_h, vxc_gllb_s, &
536 vxc_tmp_h, vxc_tmp_s, vtau, dummy, e_uniform, drho_h, drho_s, vxg, atom_list, &
537 atomic_kind_set, qs_kind_set, deriv, atomic_grid, rho_struct_ao, &
538 harmonics, molecular_orbitals, rho_structure, r_h, r_s, dr_h, dr_s, &
539 r_h_d, r_s_d, rho_atom_set, rho_atom, para_env, &
540 mo_eigenvalues, local_rho_set, matrix_ks, &
541 orbital_density_matrix, vxc_saop_h, vxc_saop_s)
542
543 ! tau is needed for fill_rho_set, but should never be used
544 NULLIFY (tau)
545 NULLIFY (dft_control, oce, sab)
546
547 CALL get_qs_env(qs_env, input=input, &
548 rho=rho_structure, &
549 mos=molecular_orbitals, &
550 atomic_kind_set=atomic_kind_set, &
551 qs_kind_set=qs_kind_set, &
552 rho_atom_set=rho_atom_set, &
553 matrix_ks=matrix_ks, &
554 dft_control=dft_control, &
555 para_env=para_env, &
556 oce=oce, sab_orb=sab)
557
558 CALL qs_rho_get(rho_structure, rho_ao=rho_struct_ao)
559
560 xc_section_orig => section_vals_get_subs_vals(input, "DFT%XC")
561 CALL section_vals_retain(xc_section_orig)
562 CALL section_vals_duplicate(xc_section_orig, xc_section_tmp)
563
564 accint = dft_control%qs_control%gapw_control%accurate_xcint
565
566 ! [SC] the following code can be traced back to SVN rev. 4296 (git:f97138b) that
567 ! has removed the component 'nspins' from the derived type 'dft_control_type'.
568 ! Is it worth to remove the code below in favour of 'dft_control%nspins'
569 ! since its reintroduction? Note that in case of ROKS calculations,
570 ! 'lsd == .FALSE.' but 'dft_control%nspins == 2'.
571 CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
572 IF (lsd) THEN
573 nspins = 2
574 ELSE
575 nspins = 1
576 END IF
577
578 CALL section_vals_val_get(xc_section_orig, "DENSITY_CUTOFF", &
579 r_val=density_cut)
580 CALL section_vals_val_get(xc_section_orig, "GRADIENT_CUTOFF", &
581 r_val=gradient_cut)
582 CALL section_vals_val_get(xc_section_orig, "TAU_CUTOFF", &
583 r_val=tau_cut)
584
585 ! remap pointer
586 ns = SIZE(rho_struct_ao)
587 psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
588 CALL calculate_rho_atom_coeff(qs_env, psmat, rho_atom_set, qs_kind_set, oce, sab, para_env)
589 CALL prepare_gapw_den(qs_env)
590
591 ALLOCATE (mo_coeff(nspins), single_mo_coeff(nspins), mo_eigenvalues(nspins))
592
593 CALL dbcsr_allocate_matrix_set(orbital_density_matrix, nspins)
594
595 DO ispin = 1, nspins
596 CALL get_mo_set(molecular_orbitals(ispin), &
597 mo_coeff=mo_coeff(ispin)%matrix, &
598 eigenvalues=mo_eigenvalues(ispin)%array, &
599 homo=homo(ispin))
600 CALL cp_fm_create(single_mo_coeff(ispin), &
601 mo_coeff(ispin)%matrix%matrix_struct, &
602 "orbital density matrix")
603 CALL cp_fm_get_info(single_mo_coeff(ispin), &
604 nrow_global=nrow(ispin), ncol_global=ncol(ispin))
605 ALLOCATE (orbital_density_matrix(ispin)%matrix)
606 CALL dbcsr_copy(orbital_density_matrix(ispin)%matrix, &
607 rho_struct_ao(ispin)%matrix, &
608 "orbital density")
609 END DO
610 CALL local_rho_set_create(local_rho_set)
611 CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
612 qs_kind_set, dft_control, para_env)
613
614 DO ikind = 1, SIZE(atomic_kind_set)
615 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
616
617 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
618 harmonics=harmonics, grid_atom=atomic_grid)
619 IF (.NOT. paw_atom) cycle
620
621 nr = atomic_grid%nr
622 na = atomic_grid%ng_sphere
623 bounds(1:2, 1:3) = 1
624 bounds(2, 1) = na
625 bounds(2, 2) = nr
626
627 CALL xc_dset_create(deriv_set, local_bounds=bounds)
628
629 CALL xc_rho_set_create(rho_set_h, bounds, density_cut, &
630 gradient_cut, tau_cut)
631 CALL xc_rho_set_create(rho_set_s, bounds, density_cut, &
632 gradient_cut, tau_cut)
633 CALL xc_rho_set_create(orb_rho_set_h, bounds, density_cut, &
634 gradient_cut, tau_cut)
635 CALL xc_rho_set_create(orb_rho_set_s, bounds, density_cut, &
636 gradient_cut, tau_cut)
637
638 CALL xc_rho_cflags_setall(needs, .false.)
639 IF (lsd) THEN
640 CALL xb88_lsd_info(needs=needs)
641 needs%norm_drho = .true.
642 ELSE
643 CALL xb88_lda_info(needs=needs)
644 END IF
645 CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
646 CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
647 CALL xc_rho_cflags_setall(needs_orbs, .false.)
648 needs_orbs%rho = .true.
649 IF (lsd) needs_orbs%rho_spin = .true.
650 CALL xc_rho_set_atom_update(orb_rho_set_h, needs, nspins, bounds)
651 CALL xc_rho_set_atom_update(orb_rho_set_s, needs, nspins, bounds)
652
653 ALLOCATE (rho_h(1:na, 1:nr, 1:nspins), rho_s(1:na, 1:nr, 1:nspins))
654 ALLOCATE (weight_h(1:na, 1:nr), weight_s(1:na, 1:nr))
655 ALLOCATE (vxc_lb_h(1:na, 1:nr, 1:nspins), vxc_lb_s(1:na, 1:nr, 1:nspins))
656 ALLOCATE (vxc_gllb_h(1:na, 1:nr, 1:nspins), vxc_gllb_s(1:na, 1:nr, 1:nspins))
657 ALLOCATE (vxc_tmp_h(1:na, 1:nr, 1:nspins), vxc_tmp_s(1:na, 1:nr, 1:nspins))
658 ALLOCATE (vxc_saop_h(1:na, 1:nr, 1:nspins), vxc_saop_s(1:na, 1:nr, 1:nspins))
659 ALLOCATE (drho_h(1:4, 1:na, 1:nr, 1:nspins), drho_s(1:4, 1:na, 1:nr, 1:nspins))
660
661 ! Distribute the atoms of this kind
662 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
663
664 DO ir = 1, nr
665 DO ia = 1, na
666 weight_h(ia, ir) = atomic_grid%wr(ir)*atomic_grid%wa(ia)
667 END DO
668 END DO
669 IF (accint) THEN
670 on = dft_control%qs_control%gapw_control%oweights
671 alpha = dft_control%qs_control%gapw_control%aw(ikind)
672 ALLOCATE (fw(nr))
673 CALL calc_weight_function(fw, atomic_grid%rad2, on, alpha)
674 DO ir = 1, nr
675 agr = 1.0_dp - fw(ir)
676 DO ia = 1, na
677 weight_s(ia, ir) = agr*atomic_grid%wr(ir)*atomic_grid%wa(ia)
678 END DO
679 END DO
680 DEALLOCATE (fw)
681 ELSE
682 weight_s(:, :) = weight_h(:, :)
683 END IF
684
685 DO iat = 1, natom !bo(1),bo(2)
686 iatom = atom_list(iat)
687
688 rho_atom => rho_atom_set(iatom)
689 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
690 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
691 rho_rad_s=r_s, drho_rad_h=dr_h, &
692 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
693 rho_rad_s_d=r_s_d)
694 rho_h = 0.0_dp
695 rho_s = 0.0_dp
696 drho_h = 0.0_dp
697 drho_s = 0.0_dp
698 DO ir = 1, nr
699 CALL calc_rho_angular(atomic_grid, harmonics, nspins, .true., &
700 ir, r_h, r_s, rho_h, rho_s, &
701 dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
702 END DO
703 DO ir = 1, nr
704 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau, na, ir)
705 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau, na, ir)
706 END DO
707
708 !-----------------------------!
709 ! 1. Slater exchange for LB94 !
710 !-----------------------------!
711 xc_fun_section_orig => section_vals_get_subs_vals(xc_section_orig, &
712 "XC_FUNCTIONAL")
713 CALL section_vals_create(xc_fun_section_tmp, xc_fun_section_orig%section)
714 CALL section_vals_val_set(xc_fun_section_tmp, "_SECTION_PARAMETERS_", &
715 i_val=xc_funct_no_shortcut)
716 CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
717 l_val=.true.)
718 CALL section_vals_set_subs_vals(xc_section_tmp, "XC_FUNCTIONAL", &
719 xc_fun_section_tmp)
720
721 !---------------------!
722 ! Both: hard and soft !
723 !---------------------!
724 CALL xc_dset_zero_all(deriv_set)
725 CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_h, deriv_set, 1, needs, &
726 weight_h, lsd, na, nr, exc, vxc_tmp_h, vxg, vtau)
727 CALL xc_dset_zero_all(deriv_set)
728 CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_s, deriv_set, 1, needs, &
729 weight_s, lsd, na, nr, exc, vxc_tmp_s, vxg, vtau)
730
731 !-------------------------------------------!
732 ! 2. PZ correlation for LB94 and ec_uniform !
733 !-------------------------------------------!
734 CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
735 l_val=.false.)
736 CALL section_vals_val_set(xc_fun_section_tmp, "PZ81%_SECTION_PARAMETERS_", &
737 l_val=.true.)
738
739 !------!
740 ! Hard !
741 !------!
742 CALL xc_dset_zero_all(deriv_set)
743 CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_h, deriv_set, 1, needs, &
744 weight_h, lsd, na, nr, exc, vxc_lb_h, vxg, vtau)
745 vxc_lb_h = vxc_lb_h + alpha*vxc_tmp_h
746 DO ispin = 1, nspins
747 dummy => vxc_tmp_h(:, :, ispin:ispin)
748 CALL add_lb_pot(dummy, rho_set_h, lsd, ispin)
749 vxc_lb_h(:, :, ispin) = vxc_lb_h(:, :, ispin) - weight_h(:, :)*vxc_tmp_h(:, :, ispin)
750 END DO
751 NULLIFY (dummy)
752
753 vxc_gllb_h = 0.0_dp
754 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
755 cpassert(ASSOCIATED(deriv))
756 CALL xc_derivative_get(deriv, deriv_data=e_uniform)
757 DO ispin = 1, nspins
758 dummy => vxc_gllb_h(:, :, ispin:ispin)
759 CALL calc_2excpbe(dummy, rho_set_h, e_uniform, lsd)
760 vxc_gllb_h(:, :, ispin) = vxc_gllb_h(:, :, ispin)*weight_h(:, :)
761 END DO
762 NULLIFY (deriv, dummy, e_uniform)
763
764 !------!
765 ! Soft !
766 !------!
767 CALL xc_dset_zero_all(deriv_set)
768 CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_s, deriv_set, 1, needs, &
769 weight_s, lsd, na, nr, exc, vxc_lb_s, vxg, vtau)
770
771 vxc_lb_s = vxc_lb_s + alpha*vxc_tmp_s
772 DO ispin = 1, nspins
773 dummy => vxc_tmp_s(:, :, ispin:ispin)
774 CALL add_lb_pot(dummy, rho_set_s, lsd, ispin)
775 vxc_lb_s(:, :, ispin) = vxc_lb_s(:, :, ispin) - weight_s(:, :)*vxc_tmp_s(:, :, ispin)
776 END DO
777 NULLIFY (dummy)
778
779 vxc_gllb_s = 0.0_dp
780 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
781 cpassert(ASSOCIATED(deriv))
782 CALL xc_derivative_get(deriv, deriv_data=e_uniform)
783 DO ispin = 1, nspins
784 dummy => vxc_gllb_s(:, :, ispin:ispin)
785 CALL calc_2excpbe(dummy, rho_set_s, e_uniform, lsd)
786 vxc_gllb_s(:, :, ispin) = vxc_gllb_s(:, :, ispin)*weight_s(:, :)
787 END DO
788 NULLIFY (deriv, dummy, e_uniform)
789
790 !------------------!
791 ! Now the orbitals !
792 !------------------!
793 vxc_tmp_h = 0.0_dp; vxc_tmp_s = 0.0_dp
794
795 DO ispin = 1, nspins
796
797 DO orb = 1, homo(ispin) - 1
798
799 ALLOCATE (coeff_col(nrow(ispin), 1))
800
801 efac = k_rho*sqrt(mo_eigenvalues(ispin)%array(homo(ispin)) - &
802 mo_eigenvalues(ispin)%array(orb))
803 IF (.NOT. lsd) efac = 2.0_dp*efac
804
805 CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
806 CALL cp_fm_get_submatrix(mo_coeff(ispin)%matrix, coeff_col, &
807 1, orb, nrow(ispin), 1)
808 CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
809 1, orb)
810 CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
811 CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
812 matrix_v=single_mo_coeff(ispin), &
813 ncol=ncol(ispin), &
814 alpha=1.0_dp)
815
816 DEALLOCATE (coeff_col)
817
818 ! This calculates the CPC and density on the grids for every atom even though
819 ! we need it only for iatom at the moment. It seems that to circumvent this,
820 ! the routines must be adapted to calculate just iatom
821 ! remap pointer
822 ns = SIZE(orbital_density_matrix)
823 psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
824 CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
825 CALL prepare_gapw_den(qs_env, local_rho_set, .false.)
826
827 rho_atom => local_rho_set%rho_atom_set(iatom)
828 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
829 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
830 rho_h = 0.0_dp
831 rho_s = 0.0_dp
832 drho_h = 0.0_dp
833 drho_s = 0.0_dp
834 DO ir = 1, nr
835 CALL calc_rho_angular(atomic_grid, harmonics, nspins, .false., &
836 ir, r_h, r_s, rho_h, rho_s, &
837 dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
838 END DO
839 DO ir = 1, nr
840 CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
841 CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
842 END DO
843
844 IF (lsd) THEN
845 IF (ispin == 1) THEN
846 vxc_tmp_h(:, :, 1) = vxc_tmp_h(:, :, 1) + efac*orb_rho_set_h%rhoa(:, :, 1)
847 vxc_tmp_s(:, :, 1) = vxc_tmp_s(:, :, 1) + efac*orb_rho_set_s%rhoa(:, :, 1)
848 ELSE
849 vxc_tmp_h(:, :, 2) = vxc_tmp_h(:, :, 2) + efac*orb_rho_set_h%rhob(:, :, 1)
850 vxc_tmp_s(:, :, 2) = vxc_tmp_s(:, :, 2) + efac*orb_rho_set_s%rhob(:, :, 1)
851 END IF
852 ELSE
853 vxc_tmp_h(:, :, 1) = vxc_tmp_h(:, :, 1) + efac*orb_rho_set_h%rho(:, :, 1)
854 vxc_tmp_s(:, :, 1) = vxc_tmp_s(:, :, 1) + efac*orb_rho_set_s%rho(:, :, 1)
855 END IF
856
857 END DO ! orb
858
859 END DO ! ispin
860
861 IF (lsd) THEN
862 DO ir = 1, nr
863 DO ia = 1, na
864 IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
865 vxc_gllb_h(ia, ir, 1) = vxc_gllb_h(ia, ir, 1) + &
866 weight_h(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
867 END IF
868 IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
869 vxc_gllb_h(ia, ir, 2) = vxc_gllb_h(ia, ir, 2) + &
870 weight_h(ia, ir)*vxc_tmp_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
871 END IF
872 IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
873 vxc_gllb_s(ia, ir, 1) = vxc_gllb_s(ia, ir, 1) + &
874 weight_s(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
875 END IF
876 IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
877 vxc_gllb_s(ia, ir, 2) = vxc_gllb_s(ia, ir, 2) + &
878 weight_s(ia, ir)*vxc_tmp_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
879 END IF
880 END DO
881 END DO
882 ELSE
883 DO ir = 1, nr
884 DO ia = 1, na
885 IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
886 vxc_gllb_h(ia, ir, 1) = vxc_gllb_h(ia, ir, 1) + &
887 weight_h(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
888 END IF
889 IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
890 vxc_gllb_s(ia, ir, 1) = vxc_gllb_s(ia, ir, 1) + &
891 weight_s(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
892 END IF
893 END DO
894 END DO
895 END IF
896
897 vxc_saop_h = 0.0_dp; vxc_saop_s = 0.0_dp
898
899 DO ispin = 1, nspins
900
901 DO orb = 1, homo(ispin)
902
903 ALLOCATE (coeff_col(nrow(ispin), 1))
904
905 we_lb = exp(-2.0_dp*(mo_eigenvalues(ispin)%array(homo(ispin)) - &
906 mo_eigenvalues(ispin)%array(orb))**2)
907 we_gllb = 1.0_dp - we_lb
908 IF (.NOT. lsd) THEN
909 we_lb = 2.0_dp*we_lb
910 we_gllb = 2.0_dp*we_gllb
911 END IF
912
913 vxc_tmp_h(:, :, ispin) = we_lb*vxc_lb_h(:, :, ispin) + &
914 we_gllb*vxc_gllb_h(:, :, ispin)
915 vxc_tmp_s(:, :, ispin) = we_lb*vxc_lb_s(:, :, ispin) + &
916 we_gllb*vxc_gllb_s(:, :, ispin)
917
918 CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
919 CALL cp_fm_get_submatrix(mo_coeff(ispin)%matrix, coeff_col, &
920 1, orb, nrow(ispin), 1)
921 CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
922 1, orb)
923 CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
924 CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
925 matrix_v=single_mo_coeff(ispin), &
926 ncol=ncol(ispin), &
927 alpha=1.0_dp)
928
929 DEALLOCATE (coeff_col)
930
931 ! This calculates the CPC and density on the grids for every atom even though
932 ! we need it only for iatom at the moment. It seems that to circumvent this,
933 ! the routines must be adapted to calculate just iatom
934 ! remap pointer
935 ns = SIZE(orbital_density_matrix)
936 psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
937 CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
938 CALL prepare_gapw_den(qs_env, local_rho_set, .false.)
939
940 rho_atom => local_rho_set%rho_atom_set(iatom)
941 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
942 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
943 rho_h = 0.0_dp
944 rho_s = 0.0_dp
945 drho_h = 0.0_dp
946 drho_s = 0.0_dp
947 DO ir = 1, nr
948 CALL calc_rho_angular(atomic_grid, harmonics, nspins, .false., &
949 ir, r_h, r_s, rho_h, rho_s, &
950 dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
951 END DO
952 DO ir = 1, nr
953 CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
954 CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
955 END DO
956
957 IF (lsd) THEN
958 IF (ispin == 1) THEN
959 vxc_saop_h(:, :, 1) = vxc_saop_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rhoa(:, :, 1)
960 vxc_saop_s(:, :, 1) = vxc_saop_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rhoa(:, :, 1)
961 ELSE
962 vxc_saop_h(:, :, 2) = vxc_saop_h(:, :, 2) + vxc_tmp_h(:, :, 2)*orb_rho_set_h%rhob(:, :, 1)
963 vxc_saop_s(:, :, 2) = vxc_saop_s(:, :, 2) + vxc_tmp_s(:, :, 2)*orb_rho_set_s%rhob(:, :, 1)
964 END IF
965 ELSE
966 vxc_saop_h(:, :, 1) = vxc_saop_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rho(:, :, 1)
967 vxc_saop_s(:, :, 1) = vxc_saop_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rho(:, :, 1)
968 END IF
969
970 END DO ! orb
971
972 END DO ! ispin
973
974 IF (lsd) THEN
975 DO ir = 1, nr
976 DO ia = 1, na
977 IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
978 vxc_saop_h(ia, ir, 1) = vxc_saop_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
979 ELSE
980 vxc_saop_h(ia, ir, 1) = 0.0_dp
981 END IF
982 IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
983 vxc_saop_h(ia, ir, 2) = vxc_saop_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
984 ELSE
985 vxc_saop_h(ia, ir, 2) = 0.0_dp
986 END IF
987 IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
988 vxc_saop_s(ia, ir, 1) = vxc_saop_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
989 ELSE
990 vxc_saop_s(ia, ir, 1) = 0.0_dp
991 END IF
992 IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
993 vxc_saop_s(ia, ir, 2) = vxc_saop_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
994 ELSE
995 vxc_saop_s(ia, ir, 2) = 0.0_dp
996 END IF
997 END DO
998 END DO
999 ELSE
1000 DO ir = 1, nr
1001 DO ia = 1, na
1002 IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
1003 vxc_saop_h(ia, ir, 1) = vxc_saop_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
1004 ELSE
1005 vxc_saop_h(ia, ir, 1) = 0.0_dp
1006 END IF
1007 IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
1008 vxc_saop_s(ia, ir, 1) = vxc_saop_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
1009 ELSE
1010 vxc_saop_s(ia, ir, 1) = 0.0_dp
1011 END IF
1012 END DO
1013 END DO
1014 END IF
1015
1016 rho_atom => rho_atom_set(iatom)
1017 CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1018 CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis, &
1019 harmonics=harmonics, grid_atom=grid_atom)
1020 SELECT CASE (oe_corr)
1021 CASE (oe_lb)
1022 CALL gavxcgb_nogc(vxc_lb_h, vxc_lb_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1023 CASE (oe_gllb)
1024 CALL gavxcgb_nogc(vxc_gllb_h, vxc_gllb_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1025 CASE (oe_saop)
1026 CALL gavxcgb_nogc(vxc_saop_h, vxc_saop_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1027 CASE default
1028 cpabort("Unknown correction!")
1029 END SELECT
1030
1031 END DO
1032
1033 DEALLOCATE (rho_h, rho_s, weight_h, weight_s)
1034 DEALLOCATE (vxc_lb_h, vxc_lb_s)
1035 DEALLOCATE (vxc_gllb_h, vxc_gllb_s)
1036 DEALLOCATE (vxc_tmp_h, vxc_tmp_s)
1037 DEALLOCATE (vxc_saop_h, vxc_saop_s)
1038 DEALLOCATE (drho_h, drho_s)
1039
1040 CALL xc_dset_release(deriv_set)
1041 CALL xc_rho_set_release(rho_set_h)
1042 CALL xc_rho_set_release(rho_set_s)
1043 CALL xc_rho_set_release(orb_rho_set_h)
1044 CALL xc_rho_set_release(orb_rho_set_s)
1045
1046 END DO
1047
1048 ! remap pointer
1049 ns = SIZE(matrix_ks)
1050 ksmat(1:ns, 1:1) => matrix_ks(1:ns)
1051 ns = SIZE(rho_struct_ao)
1052 psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
1053
1054 CALL update_ks_atom(qs_env, ksmat, psmat, forces=.false.)
1055
1056 !---------!
1057 ! Cleanup !
1058 !---------!
1059 CALL section_vals_release(xc_fun_section_tmp)
1060 CALL section_vals_release(xc_section_tmp)
1061 CALL section_vals_release(xc_section_orig)
1062
1063 CALL local_rho_set_release(local_rho_set)
1064 CALL cp_fm_release(single_mo_coeff)
1065 DEALLOCATE (mo_coeff, mo_eigenvalues)
1066 CALL dbcsr_deallocate_matrix_set(orbital_density_matrix)
1067
1068 END SUBROUTINE gapw_add_atomic_saop_pot
1069
1070! **************************************************************************************************
1071!> \brief ...
1072!> \param pot ...
1073!> \param rho_set ...
1074!> \param lsd ...
1075!> \param spin ...
1076! **************************************************************************************************
1077 SUBROUTINE add_lb_pot(pot, rho_set, lsd, spin)
1078
1079 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pot
1080 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1081 LOGICAL, INTENT(IN) :: lsd
1082 INTEGER, INTENT(IN) :: spin
1083
1084 REAL(kind=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp
1085
1086 INTEGER :: i, j, k
1087 INTEGER, DIMENSION(2, 3) :: bo
1088 REAL(kind=dp) :: n, n_13, x, x2
1089
1090 bo = rho_set%local_bounds
1091
1092 DO k = bo(1, 3), bo(2, 3)
1093 DO j = bo(1, 2), bo(2, 2)
1094 DO i = bo(1, 1), bo(2, 1)
1095 IF (.NOT. lsd) THEN
1096 IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1097 n = rho_set%rho(i, j, k)/2.0_dp
1098 n_13 = n**ob3
1099 x = (rho_set%norm_drho(i, j, k)/2.0_dp)/(n*n_13)
1100 x2 = x*x
1101 pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(x + sqrt(x2 + 1.0_dp)))
1102 END IF
1103 ELSE
1104 IF (spin == 1) THEN
1105 IF (rho_set%rhoa(i, j, k) > rho_set%rho_cutoff) THEN
1106 n_13 = rho_set%rhoa_1_3(i, j, k)
1107 x = rho_set%norm_drhoa(i, j, k)/(rho_set%rhoa(i, j, k)*n_13)
1108 x2 = x*x
1109 pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(sqrt(x2 + 1.0_dp) + x))
1110 END IF
1111 ELSE IF (spin == 2) THEN
1112 IF (rho_set%rhob(i, j, k) > rho_set%rho_cutoff) THEN
1113 n_13 = rho_set%rhob_1_3(i, j, k)
1114 x = rho_set%norm_drhob(i, j, k)/(rho_set%rhob(i, j, k)*n_13)
1115 x2 = x*x
1116 pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(sqrt(x2 + 1.0_dp) + x))
1117 END IF
1118 END IF
1119 END IF
1120 END DO
1121 END DO
1122 END DO
1123
1124 END SUBROUTINE add_lb_pot
1125
1126! **************************************************************************************************
1127!> \brief ...
1128!> \param pot ...
1129!> \param rho_set ...
1130!> \param e_uniform ...
1131!> \param lsd ...
1132! **************************************************************************************************
1133 SUBROUTINE calc_2excpbe(pot, rho_set, e_uniform, lsd)
1134
1135 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pot
1136 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1137 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: e_uniform
1138 LOGICAL, INTENT(IN) :: lsd
1139
1140 INTEGER :: i, j, k
1141 INTEGER, DIMENSION(2, 3) :: bo
1142 REAL(kind=dp) :: e_unif, rho
1143
1144 bo = rho_set%local_bounds
1145
1146 DO k = bo(1, 3), bo(2, 3)
1147 DO j = bo(1, 2), bo(2, 2)
1148 DO i = bo(1, 1), bo(2, 1)
1149 IF (.NOT. lsd) THEN
1150 IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1151 e_unif = e_uniform(i, j, k)/rho_set%rho(i, j, k)
1152 ELSE
1153 e_unif = 0.0_dp
1154 END IF
1155 pot(i, j, k) = &
1156 2.0_dp* &
1157 calc_ecpbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1158 e_unif, rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1159 2.0_dp* &
1160 calc_expbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1161 rho_set%rho_cutoff, rho_set%drho_cutoff)
1162 ELSE
1163 rho = rho_set%rhoa(i, j, k) + rho_set%rhob(i, j, k)
1164 IF (rho > rho_set%rho_cutoff) THEN
1165 e_unif = e_uniform(i, j, k)/rho
1166 ELSE
1167 e_unif = 0.0_dp
1168 END IF
1169 pot(i, j, k) = &
1170 2.0_dp* &
1171 calc_ecpbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1172 e_unif, &
1173 rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1174 2.0_dp* &
1175 calc_expbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1176 rho_set%rho_cutoff, rho_set%drho_cutoff)
1177 END IF
1178 END DO
1179 END DO
1180 END DO
1181
1182 END SUBROUTINE calc_2excpbe
1183
1184! **************************************************************************************************
1185!> \brief ...
1186!> \param ra ...
1187!> \param rb ...
1188!> \param ngr ...
1189!> \param ec_unif ...
1190!> \param rc ...
1191!> \param ngrc ...
1192!> \return ...
1193! **************************************************************************************************
1194 FUNCTION calc_ecpbe_u(ra, rb, ngr, ec_unif, rc, ngrc) RESULT(res)
1195
1196 REAL(kind=dp), INTENT(in) :: ra, rb, ngr, ec_unif, rc, ngrc
1197 REAL(kind=dp) :: res
1198
1199 REAL(kind=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp, tb3 = 2.0_dp/3.0_dp
1200
1201 REAL(kind=dp) :: a, at2, h, kf, kl, ks, phi, phi3, r, t2, &
1202 zeta
1203
1204 r = ra + rb
1205 h = 0.0_dp
1206 IF (r > rc .AND. ngr > ngrc) THEN
1207 zeta = (ra - rb)/r
1208 IF (zeta > 1.0_dp) zeta = 1.0_dp ! machine precision problem
1209 IF (zeta < -1.0_dp) zeta = -1.0_dp ! machine precision problem
1210 phi = ((1.0_dp + zeta)**tb3 + (1.0_dp - zeta)**tb3)/2.0_dp
1211 phi3 = phi*phi*phi
1212 kf = (3.0_dp*r*pi*pi)**ob3
1213 ks = sqrt(4.0_dp*kf/pi)
1214 t2 = (ngr/(2.0_dp*phi*ks*r))**2
1215 a = beta_ec/gamma_saop/(exp(-ec_unif/(gamma_saop*phi3)) - 1.0_dp)
1216 at2 = a*t2
1217 kl = (1.0_dp + at2)/(1.0_dp + at2 + at2*at2)
1218 h = gamma_saop*log(1.0_dp + beta_ec/gamma_saop*t2*kl)
1219 END IF
1220 res = ec_unif + h
1221
1222 END FUNCTION calc_ecpbe_u
1223
1224! **************************************************************************************************
1225!> \brief ...
1226!> \param r ...
1227!> \param ngr ...
1228!> \param ec_unif ...
1229!> \param rc ...
1230!> \param ngrc ...
1231!> \return ...
1232! **************************************************************************************************
1233 FUNCTION calc_ecpbe_r(r, ngr, ec_unif, rc, ngrc) RESULT(res)
1234
1235 REAL(kind=dp), INTENT(in) :: r, ngr, ec_unif, rc, ngrc
1236 REAL(kind=dp) :: res
1237
1238 REAL(kind=dp) :: a, at2, h, kf, kl, ks, t2
1239
1240 h = 0.0_dp
1241 IF (r > rc .AND. ngr > ngrc) THEN
1242 kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1243 ks = sqrt(4.0_dp*kf/pi)
1244 t2 = (ngr/(2.0_dp*ks*r))**2
1245 a = beta_ec/gamma_saop/(exp(-ec_unif/gamma_saop) - 1.0_dp)
1246 at2 = a*t2
1247 kl = (1.0_dp + at2)/(1.0_dp + at2 + at2*at2)
1248 h = gamma_saop*log(1.0_dp + beta_ec/gamma_saop*t2*kl)
1249 END IF
1250 res = ec_unif + h
1251
1252 END FUNCTION calc_ecpbe_r
1253
1254! **************************************************************************************************
1255!> \brief ...
1256!> \param ra ...
1257!> \param rb ...
1258!> \param ngr ...
1259!> \param rc ...
1260!> \param ngrc ...
1261!> \return ...
1262! **************************************************************************************************
1263 FUNCTION calc_expbe_u(ra, rb, ngr, rc, ngrc) RESULT(res)
1264
1265 REAL(kind=dp), INTENT(in) :: ra, rb, ngr, rc, ngrc
1266 REAL(kind=dp) :: res
1267
1268 REAL(kind=dp) :: r
1269
1270 r = ra + rb
1271 res = calc_expbe_r(r, ngr, rc, ngrc)
1272
1273 END FUNCTION calc_expbe_u
1274
1275! **************************************************************************************************
1276!> \brief ...
1277!> \param r ...
1278!> \param ngr ...
1279!> \param rc ...
1280!> \param ngrc ...
1281!> \return ...
1282! **************************************************************************************************
1283 FUNCTION calc_expbe_r(r, ngr, rc, ngrc) RESULT(res)
1284
1285 REAL(kind=dp), INTENT(in) :: r, ngr, rc, ngrc
1286 REAL(kind=dp) :: res
1287
1288 REAL(kind=dp) :: ex_unif, fx, kf, s
1289
1290 IF (r > rc) THEN
1291 kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1292 ex_unif = -3.0_dp*kf/(4.0_dp*pi)
1293 fx = 1.0_dp
1294 IF (ngr > ngrc) THEN
1295 s = ngr/(2.0_dp*kf*r)
1296 fx = fx + kappa - kappa/(1.0_dp + mu*s*s/kappa)
1297 END IF
1298 res = ex_unif*fx
1299 ELSE
1300 res = 0.0_dp
1301 END IF
1302
1303 END FUNCTION calc_expbe_r
1304
1305END MODULE xc_pot_saop
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:95
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:103
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:89
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:974
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:1045
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:757
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:535
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:166
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:932
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:629
input_constants::xc_funct_no_shortcut
integer, parameter, public xc_funct_no_shortcut
Definition input_constants.F:593
input_constants::do_method_gapw
integer, parameter, public do_method_gapw
Definition input_constants.F:270
input_constants::oe_lb
integer, parameter, public oe_lb
Definition input_constants.F:629
input_constants::oe_gllb
integer, parameter, public oe_gllb
Definition input_constants.F:629
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:27
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:1713
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, matrix_vhxc, 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:544
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, hund_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, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, 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:491
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, cmo_coeff)
Get the components of a MO set data structure.
Definition qs_mo_types.F:445
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:1023
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:560
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_utils
Support routines for integrals of the Vxc/Fxc/Gxc potentials calculated for the atomic density in the...
Definition qs_vxc_atom_utils.F:12
qs_vxc_atom_utils::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_utils.F:219
qs_vxc_atom_utils::gavxcgb_nogc
subroutine, public gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition qs_vxc_atom_utils.F:595
qs_vxc_atom_utils::calc_weight_function
subroutine, public calc_weight_function(fun, r2, order, alpha)
Calculates the radial weight function for accurate XC integration.
Definition qs_vxc_atom_utils.F:1108
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:523
xc_atom::fill_rho_set
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition xc_atom.F:683
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:739
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:127
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:478
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:786
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:129
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:115
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:768
orbital
Orbital angular momentum.
Definition grid_common.h:134
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:223
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:190
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:53
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