Line data Source code
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 : ! **************************************************************************************************
11 : MODULE 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 :
117 : CONTAINS
118 :
119 : ! **************************************************************************************************
120 : !> \brief ...
121 : !> \param ks_matrix ...
122 : !> \param qs_env ...
123 : !> \param oe_corr ...
124 : ! **************************************************************************************************
125 14 : 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 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_col
140 : REAL(KIND=dp), DIMENSION(3, 3) :: virial_xc_tmp
141 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
142 14 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: e_uniform
143 14 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: single_mo_coeff
144 : TYPE(cp_fm_type), POINTER :: mo_coeff
145 28 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: orbital_density_matrix, rho_struct_ao
146 14 : TYPE(mo_set_type), DIMENSION(:), POINTER :: molecular_orbitals
147 : TYPE(pw_c1d_gs_type) :: orbital_g
148 14 : 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 14 : TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: vxc_GLLB, vxc_SAOP
153 28 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_struct_r, tau, vxc_LB, &
154 14 : 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 14 : NULLIFY (ks_env, pw_env, auxbas_pw_pool, input)
168 14 : NULLIFY (rho_g, rho_r, tau, rho_struct, e_uniform)
169 14 : NULLIFY (vxc_LB, vxc_tmp, vxc_tau)
170 14 : NULLIFY (mo_eigenvalues, deriv, rho_struct_r, rho_struct_ao)
171 14 : 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 14 : mos=molecular_orbitals)
181 14 : compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
182 14 : CALL section_vals_val_get(input, "DFT%QS%METHOD", i_val=dft_method_id)
183 14 : gapw = (dft_method_id == do_method_gapw)
184 :
185 14 : xc_section_orig => section_vals_get_subs_vals(input, "DFT%XC")
186 14 : CALL section_vals_retain(xc_section_orig)
187 14 : CALL section_vals_duplicate(xc_section_orig, xc_section_tmp)
188 :
189 : CALL section_vals_val_get(xc_section_orig, "DENSITY_CUTOFF", &
190 14 : r_val=density_cut)
191 : CALL section_vals_val_get(xc_section_orig, "GRADIENT_CUTOFF", &
192 14 : r_val=gradient_cut)
193 : CALL section_vals_val_get(xc_section_orig, "TAU_CUTOFF", &
194 14 : r_val=tau_cut)
195 :
196 14 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
197 :
198 14 : CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
199 14 : IF (lsd) THEN
200 6 : nspins = 2
201 : ELSE
202 8 : nspins = 1
203 : END IF
204 :
205 62 : ALLOCATE (single_mo_coeff(nspins))
206 14 : CALL dbcsr_allocate_matrix_set(orbital_density_matrix, nspins)
207 14 : CALL qs_rho_get(rho_struct, rho_r=rho_struct_r, rho_ao=rho_struct_ao)
208 14 : rho_r => rho_struct_r
209 34 : DO ispin = 1, nspins
210 20 : ALLOCATE (orbital_density_matrix(ispin)%matrix)
211 : CALL dbcsr_copy(orbital_density_matrix(ispin)%matrix, &
212 34 : rho_struct_ao(ispin)%matrix, "orbital density")
213 : END DO
214 140 : 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 14 : tau_cut)
223 14 : CALL xc_rho_cflags_setall(needs, .FALSE.)
224 14 : IF (lsd) THEN
225 6 : CALL xb88_lsd_info(needs=needs)
226 6 : needs%norm_drho = .TRUE.
227 : ELSE
228 8 : CALL xb88_lda_info(needs=needs)
229 : END IF
230 : CALL section_vals_val_get(xc_section_orig, "XC_GRID%XC_DERIV", &
231 14 : i_val=xc_deriv_method_id)
232 : CALL section_vals_val_get(xc_section_orig, "XC_GRID%XC_SMOOTH_RHO", &
233 14 : 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 14 : 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 14 : "XC_FUNCTIONAL")
244 14 : 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 14 : i_val=xc_funct_no_shortcut)
247 : CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
248 14 : l_val=.TRUE.)
249 : CALL section_vals_set_subs_vals(xc_section_tmp, "XC_FUNCTIONAL", &
250 14 : xc_fun_section_tmp)
251 :
252 14 : 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 14 : compute_virial=.FALSE., virial_xc=virial_xc_tmp)
256 :
257 : CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
258 14 : l_val=.FALSE.)
259 : CALL section_vals_val_set(xc_fun_section_tmp, "PZ81%_SECTION_PARAMETERS_", &
260 14 : l_val=.TRUE.)
261 :
262 14 : 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 14 : compute_virial=.FALSE., virial_xc=virial_xc_tmp)
266 :
267 34 : DO ispin = 1, nspins
268 34 : CALL pw_axpy(vxc_tmp(ispin), vxc_LB(ispin), alpha)
269 : END DO
270 :
271 34 : DO ispin = 1, nspins
272 20 : CALL add_lb_pot(vxc_tmp(ispin)%array, rho_set, lsd, ispin)
273 34 : 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 14 : 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 14 : deriv_order=0)
285 :
286 14 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
287 14 : CALL xc_derivative_get(deriv, deriv_data=e_uniform)
288 :
289 62 : ALLOCATE (vxc_GLLB(nspins))
290 34 : DO ispin = 1, nspins
291 34 : CALL auxbas_pw_pool%create_pw(vxc_GLLB(ispin))
292 : END DO
293 :
294 34 : DO ispin = 1, nspins
295 34 : CALL calc_2excpbe(vxc_GLLB(ispin)%array, rho_set, e_uniform, lsd)
296 : END DO
297 :
298 14 : CALL xc_dset_release(deriv_set)
299 :
300 14 : CALL auxbas_pw_pool%create_pw(orbital)
301 14 : CALL auxbas_pw_pool%create_pw(orbital_g)
302 :
303 34 : DO ispin = 1, nspins
304 :
305 : CALL get_mo_set(molecular_orbitals(ispin), &
306 : mo_coeff=mo_coeff, &
307 : eigenvalues=mo_eigenvalues, &
308 20 : homo=homo)
309 : CALL cp_fm_create(single_mo_coeff(ispin), &
310 : mo_coeff%matrix_struct, &
311 20 : "orbital density matrix")
312 :
313 : CALL cp_fm_get_info(single_mo_coeff(ispin), &
314 20 : nrow_global=nrow(ispin), ncol_global=ncol(ispin))
315 60 : ALLOCATE (coeff_col(nrow(ispin), 1))
316 :
317 20 : CALL pw_zero(vxc_tmp(ispin))
318 :
319 98 : DO orb = 1, homo - 1
320 :
321 78 : efac = K_rho*SQRT(mo_eigenvalues(homo) - mo_eigenvalues(orb))
322 78 : IF (.NOT. lsd) efac = 2.0_dp*efac
323 :
324 78 : CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
325 : CALL cp_fm_get_submatrix(mo_coeff, coeff_col, &
326 78 : 1, orb, nrow(ispin), 1)
327 : CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
328 78 : 1, orb)
329 78 : 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 78 : alpha=1.0_dp)
334 78 : CALL pw_zero(orbital)
335 78 : 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 78 : ks_env=ks_env)
339 :
340 98 : CALL pw_axpy(orbital, vxc_tmp(ispin), efac)
341 :
342 : END DO
343 20 : DEALLOCATE (coeff_col)
344 :
345 694 : DO k = bo(1, 3), bo(2, 3)
346 24228 : DO j = bo(1, 2), bo(2, 2)
347 447395 : DO i = bo(1, 1), bo(2, 1)
348 446721 : 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 423187 : rho_r(ispin)%array(i, j, k)
351 : ELSE
352 0 : vxc_tmp(ispin)%array(i, j, k) = 0.0_dp
353 : END IF
354 : END DO
355 : END DO
356 : END DO
357 :
358 54 : CALL pw_axpy(vxc_tmp(ispin), vxc_GLLB(ispin), 1.0_dp)
359 :
360 : END DO
361 :
362 : !---------------!
363 : ! Assemble SAOP !
364 : !---------------!
365 48 : ALLOCATE (vxc_SAOP(nspins))
366 :
367 34 : DO ispin = 1, nspins
368 :
369 : CALL get_mo_set(molecular_orbitals(ispin), &
370 : mo_coeff=mo_coeff, &
371 : eigenvalues=mo_eigenvalues, &
372 20 : homo=homo)
373 20 : CALL auxbas_pw_pool%create_pw(vxc_SAOP(ispin))
374 20 : CALL pw_zero(vxc_SAOP(ispin))
375 :
376 60 : ALLOCATE (coeff_col(nrow(ispin), 1))
377 :
378 118 : DO orb = 1, homo
379 :
380 98 : we_LB = EXP(-2.0_dp*(mo_eigenvalues(homo) - mo_eigenvalues(orb))**2)
381 98 : we_GLLB = 1.0_dp - we_LB
382 98 : IF (.NOT. lsd) THEN
383 32 : we_LB = 2.0_dp*we_LB
384 32 : we_GLLB = 2.0_dp*we_GLLB
385 : END IF
386 :
387 : vxc_tmp(ispin)%array = we_LB*vxc_LB(ispin)%array + &
388 4682538 : we_GLLB*vxc_GLLB(ispin)%array
389 :
390 98 : CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
391 : CALL cp_fm_get_submatrix(mo_coeff, coeff_col, &
392 98 : 1, orb, nrow(ispin), 1)
393 : CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
394 98 : 1, orb)
395 98 : 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 98 : alpha=1.0_dp)
400 98 : CALL pw_zero(orbital)
401 98 : 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 98 : ks_env=ks_env)
405 :
406 : vxc_SAOP(ispin)%array = vxc_SAOP(ispin)%array + &
407 2341338 : orbital%array*vxc_tmp(ispin)%array
408 :
409 : END DO
410 :
411 20 : CALL dbcsr_deallocate_matrix(orbital_density_matrix(ispin)%matrix)
412 :
413 20 : DEALLOCATE (coeff_col)
414 :
415 728 : DO k = bo(1, 3), bo(2, 3)
416 24228 : DO j = bo(1, 2), bo(2, 2)
417 447395 : DO i = bo(1, 1), bo(2, 1)
418 446721 : 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 423187 : rho_r(ispin)%array(i, j, k)
421 : ELSE
422 0 : 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 14 : CALL cp_fm_release(single_mo_coeff)
431 :
432 14 : CALL xc_rho_set_release(rho_set, auxbas_pw_pool)
433 14 : CALL auxbas_pw_pool%give_back_pw(orbital)
434 14 : CALL auxbas_pw_pool%give_back_pw(orbital_g)
435 :
436 : !--------------------!
437 : ! Do the integration !
438 : !--------------------!
439 34 : DO ispin = 1, nspins
440 :
441 20 : IF (oe_corr == oe_lb) THEN
442 0 : CALL pw_copy(vxc_LB(ispin), vxc_SAOP(ispin))
443 20 : ELSE IF (oe_corr == oe_gllb) THEN
444 0 : CALL pw_copy(vxc_GLLB(ispin), vxc_SAOP(ispin))
445 : END IF
446 20 : 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 34 : gapw=gapw)
452 :
453 : END DO
454 :
455 34 : DO ispin = 1, nspins
456 20 : CALL auxbas_pw_pool%give_back_pw(vxc_SAOP(ispin))
457 20 : CALL auxbas_pw_pool%give_back_pw(vxc_GLLB(ispin))
458 20 : CALL vxc_LB(ispin)%release()
459 34 : CALL vxc_tmp(ispin)%release()
460 : END DO
461 14 : DEALLOCATE (vxc_GLLB, vxc_LB, vxc_tmp, orbital_density_matrix)
462 :
463 14 : DEALLOCATE (vxc_SAOP)
464 :
465 14 : CALL section_vals_release(xc_fun_section_tmp)
466 14 : CALL section_vals_release(xc_section_tmp)
467 14 : CALL section_vals_release(xc_section_orig)
468 :
469 : !-----------------------!
470 : ! Call the GAPW routine !
471 : !-----------------------!
472 14 : IF (gapw) THEN
473 0 : CALL gapw_add_atomic_saop_pot(qs_env, oe_corr)
474 : END IF
475 :
476 350 : END SUBROUTINE add_saop_pot
477 :
478 : ! **************************************************************************************************
479 : !> \brief ...
480 : !> \param qs_env ...
481 : !> \param oe_corr ...
482 : ! **************************************************************************************************
483 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: atom_list
493 : LOGICAL :: accint, lsd, paw_atom
494 0 : 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 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_col, weight_h, weight_s
498 0 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: dummy, e_uniform, rho_h, rho_s, vtau, &
499 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg
501 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
502 0 : TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: mo_eigenvalues
503 0 : TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff
504 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: single_mo_coeff
505 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, orbital_density_matrix, &
506 0 : rho_struct_ao
507 0 : 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 0 : 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 0 : POINTER :: sab
517 : TYPE(oce_matrix_type), POINTER :: oce
518 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
519 : TYPE(qs_rho_type), POINTER :: rho_structure
520 0 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
521 0 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
522 0 : 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 0 : NULLIFY (rho_h, rho_s, vxc_LB_h, vxc_LB_s, vxc_GLLB_h, vxc_GLLB_s, &
534 0 : vxc_tmp_h, vxc_tmp_s, vtau, dummy, e_uniform, drho_h, drho_s, vxg, atom_list, &
535 0 : atomic_kind_set, qs_kind_set, deriv, atomic_grid, rho_struct_ao, &
536 0 : harmonics, molecular_orbitals, rho_structure, r_h, r_s, dr_h, dr_s, &
537 0 : r_h_d, r_s_d, rho_atom_set, rho_atom, para_env, &
538 0 : mo_eigenvalues, local_rho_set, matrix_ks, &
539 0 : orbital_density_matrix, vxc_SAOP_h, vxc_SAOP_s)
540 :
541 : ! tau is needed for fill_rho_set, but should never be used
542 0 : NULLIFY (tau)
543 0 : 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 0 : oce=oce, sab_orb=sab)
555 :
556 0 : CALL qs_rho_get(rho_structure, rho_ao=rho_struct_ao)
557 :
558 0 : xc_section_orig => section_vals_get_subs_vals(input, "DFT%XC")
559 0 : CALL section_vals_retain(xc_section_orig)
560 0 : CALL section_vals_duplicate(xc_section_orig, xc_section_tmp)
561 :
562 0 : 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 0 : CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
570 0 : IF (lsd) THEN
571 0 : nspins = 2
572 : ELSE
573 0 : nspins = 1
574 : END IF
575 :
576 : CALL section_vals_val_get(xc_section_orig, "DENSITY_CUTOFF", &
577 0 : r_val=density_cut)
578 : CALL section_vals_val_get(xc_section_orig, "GRADIENT_CUTOFF", &
579 0 : r_val=gradient_cut)
580 : CALL section_vals_val_get(xc_section_orig, "TAU_CUTOFF", &
581 0 : r_val=tau_cut)
582 :
583 : ! remap pointer
584 0 : ns = SIZE(rho_struct_ao)
585 0 : psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
586 0 : CALL calculate_rho_atom_coeff(qs_env, psmat, rho_atom_set, qs_kind_set, oce, sab, para_env)
587 0 : CALL prepare_gapw_den(qs_env)
588 :
589 0 : ALLOCATE (mo_coeff(nspins), single_mo_coeff(nspins), mo_eigenvalues(nspins))
590 :
591 0 : CALL dbcsr_allocate_matrix_set(orbital_density_matrix, nspins)
592 :
593 0 : 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 0 : homo=homo(ispin))
598 : CALL cp_fm_create(single_mo_coeff(ispin), &
599 : mo_coeff(ispin)%matrix%matrix_struct, &
600 0 : "orbital density matrix")
601 : CALL cp_fm_get_info(single_mo_coeff(ispin), &
602 0 : nrow_global=nrow(ispin), ncol_global=ncol(ispin))
603 0 : ALLOCATE (orbital_density_matrix(ispin)%matrix)
604 : CALL dbcsr_copy(orbital_density_matrix(ispin)%matrix, &
605 : rho_struct_ao(ispin)%matrix, &
606 0 : "orbital density")
607 : END DO
608 0 : CALL local_rho_set_create(local_rho_set)
609 : CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
610 0 : qs_kind_set, dft_control, para_env)
611 :
612 0 : DO ikind = 1, SIZE(atomic_kind_set)
613 0 : 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 0 : harmonics=harmonics, grid_atom=atomic_grid)
617 0 : IF (.NOT. paw_atom) CYCLE
618 :
619 0 : nr = atomic_grid%nr
620 0 : na = atomic_grid%ng_sphere
621 0 : bounds(1:2, 1:3) = 1
622 0 : bounds(2, 1) = na
623 0 : bounds(2, 2) = nr
624 :
625 0 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
626 :
627 : CALL xc_rho_set_create(rho_set_h, bounds, density_cut, &
628 0 : gradient_cut, tau_cut)
629 : CALL xc_rho_set_create(rho_set_s, bounds, density_cut, &
630 0 : gradient_cut, tau_cut)
631 : CALL xc_rho_set_create(orb_rho_set_h, bounds, density_cut, &
632 0 : gradient_cut, tau_cut)
633 : CALL xc_rho_set_create(orb_rho_set_s, bounds, density_cut, &
634 0 : gradient_cut, tau_cut)
635 :
636 0 : CALL xc_rho_cflags_setall(needs, .FALSE.)
637 0 : IF (lsd) THEN
638 0 : CALL xb88_lsd_info(needs=needs)
639 0 : needs%norm_drho = .TRUE.
640 : ELSE
641 0 : CALL xb88_lda_info(needs=needs)
642 : END IF
643 0 : CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
644 0 : CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
645 0 : CALL xc_rho_cflags_setall(needs_orbs, .FALSE.)
646 0 : needs_orbs%rho = .TRUE.
647 0 : IF (lsd) needs_orbs%rho_spin = .TRUE.
648 0 : CALL xc_rho_set_atom_update(orb_rho_set_h, needs, nspins, bounds)
649 0 : CALL xc_rho_set_atom_update(orb_rho_set_s, needs, nspins, bounds)
650 :
651 0 : ALLOCATE (rho_h(1:na, 1:nr, 1:nspins), rho_s(1:na, 1:nr, 1:nspins))
652 0 : ALLOCATE (weight_h(1:na, 1:nr), weight_s(1:na, 1:nr))
653 0 : ALLOCATE (vxc_LB_h(1:na, 1:nr, 1:nspins), vxc_LB_s(1:na, 1:nr, 1:nspins))
654 0 : ALLOCATE (vxc_GLLB_h(1:na, 1:nr, 1:nspins), vxc_GLLB_s(1:na, 1:nr, 1:nspins))
655 0 : ALLOCATE (vxc_tmp_h(1:na, 1:nr, 1:nspins), vxc_tmp_s(1:na, 1:nr, 1:nspins))
656 0 : ALLOCATE (vxc_SAOP_h(1:na, 1:nr, 1:nspins), vxc_SAOP_s(1:na, 1:nr, 1:nspins))
657 0 : 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 0 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
661 :
662 0 : DO ir = 1, nr
663 0 : DO ia = 1, na
664 0 : weight_h(ia, ir) = atomic_grid%wr(ir)*atomic_grid%wa(ia)
665 : END DO
666 0 : IF (accint) THEN
667 0 : alpha = dft_control%qs_control%gapw_control%aw(ikind)
668 0 : agr = 1.0_dp - EXP(-alpha*atomic_grid%rad2(ir))
669 0 : DO ia = 1, na
670 0 : weight_s(ia, ir) = atomic_grid%wr(ir)*atomic_grid%wa(ia)*agr
671 : END DO
672 : ELSE
673 0 : DO ia = 1, na
674 0 : weight_s(ia, ir) = atomic_grid%wr(ir)*atomic_grid%wa(ia)
675 : END DO
676 : END IF
677 : END DO
678 :
679 0 : DO iat = 1, natom !bo(1),bo(2)
680 0 : iatom = atom_list(iat)
681 :
682 0 : rho_atom => rho_atom_set(iatom)
683 0 : 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 0 : rho_rad_s_d=r_s_d)
688 0 : rho_h = 0.0_dp
689 0 : rho_s = 0.0_dp
690 0 : drho_h = 0.0_dp
691 0 : drho_s = 0.0_dp
692 0 : 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 0 : dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
696 : END DO
697 0 : DO ir = 1, nr
698 0 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau, na, ir)
699 0 : 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 0 : "XC_FUNCTIONAL")
707 0 : 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 0 : i_val=xc_funct_no_shortcut)
710 : CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
711 0 : l_val=.TRUE.)
712 : CALL section_vals_set_subs_vals(xc_section_tmp, "XC_FUNCTIONAL", &
713 0 : xc_fun_section_tmp)
714 :
715 : !---------------------!
716 : ! Both: hard and soft !
717 : !---------------------!
718 0 : 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 0 : weight_h, lsd, na, nr, exc, vxc_tmp_h, vxg, vtau)
721 0 : 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 0 : 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 0 : l_val=.FALSE.)
730 : CALL section_vals_val_set(xc_fun_section_tmp, "PZ81%_SECTION_PARAMETERS_", &
731 0 : l_val=.TRUE.)
732 :
733 : !------!
734 : ! Hard !
735 : !------!
736 0 : 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 0 : weight_h, lsd, na, nr, exc, vxc_LB_h, vxg, vtau)
739 0 : vxc_LB_h = vxc_LB_h + alpha*vxc_tmp_h
740 0 : DO ispin = 1, nspins
741 0 : dummy => vxc_tmp_h(:, :, ispin:ispin)
742 0 : CALL add_lb_pot(dummy, rho_set_h, lsd, ispin)
743 0 : vxc_LB_h(:, :, ispin) = vxc_LB_h(:, :, ispin) - weight_h(:, :)*vxc_tmp_h(:, :, ispin)
744 : END DO
745 0 : NULLIFY (dummy)
746 :
747 0 : vxc_GLLB_h = 0.0_dp
748 0 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
749 0 : CPASSERT(ASSOCIATED(deriv))
750 0 : CALL xc_derivative_get(deriv, deriv_data=e_uniform)
751 0 : DO ispin = 1, nspins
752 0 : dummy => vxc_GLLB_h(:, :, ispin:ispin)
753 0 : CALL calc_2excpbe(dummy, rho_set_h, e_uniform, lsd)
754 0 : vxc_GLLB_h(:, :, ispin) = vxc_GLLB_h(:, :, ispin)*weight_h(:, :)
755 : END DO
756 0 : NULLIFY (deriv, dummy, e_uniform)
757 :
758 : !------!
759 : ! Soft !
760 : !------!
761 0 : 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 0 : weight_s, lsd, na, nr, exc, vxc_LB_s, vxg, vtau)
764 :
765 0 : vxc_LB_s = vxc_LB_s + alpha*vxc_tmp_s
766 0 : DO ispin = 1, nspins
767 0 : dummy => vxc_tmp_s(:, :, ispin:ispin)
768 0 : CALL add_lb_pot(dummy, rho_set_s, lsd, ispin)
769 0 : vxc_LB_s(:, :, ispin) = vxc_LB_s(:, :, ispin) - weight_s(:, :)*vxc_tmp_s(:, :, ispin)
770 : END DO
771 0 : NULLIFY (dummy)
772 :
773 0 : vxc_GLLB_s = 0.0_dp
774 0 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
775 0 : CPASSERT(ASSOCIATED(deriv))
776 0 : CALL xc_derivative_get(deriv, deriv_data=e_uniform)
777 0 : DO ispin = 1, nspins
778 0 : dummy => vxc_GLLB_s(:, :, ispin:ispin)
779 0 : CALL calc_2excpbe(dummy, rho_set_s, e_uniform, lsd)
780 0 : vxc_GLLB_s(:, :, ispin) = vxc_GLLB_s(:, :, ispin)*weight_s(:, :)
781 : END DO
782 0 : NULLIFY (deriv, dummy, e_uniform)
783 :
784 : !------------------!
785 : ! Now the orbitals !
786 : !------------------!
787 0 : vxc_tmp_h = 0.0_dp; vxc_tmp_s = 0.0_dp
788 :
789 0 : DO ispin = 1, nspins
790 :
791 0 : DO orb = 1, homo(ispin) - 1
792 :
793 0 : ALLOCATE (coeff_col(nrow(ispin), 1))
794 :
795 : efac = K_rho*SQRT(mo_eigenvalues(ispin)%array(homo(ispin)) - &
796 0 : mo_eigenvalues(ispin)%array(orb))
797 0 : IF (.NOT. lsd) efac = 2.0_dp*efac
798 :
799 0 : 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 0 : 1, orb, nrow(ispin), 1)
802 : CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
803 0 : 1, orb)
804 0 : 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 0 : alpha=1.0_dp)
809 :
810 0 : 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 0 : ns = SIZE(orbital_density_matrix)
817 0 : psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
818 0 : CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
819 0 : CALL prepare_gapw_den(qs_env, local_rho_set, .FALSE.)
820 :
821 0 : rho_atom => local_rho_set%rho_atom_set(iatom)
822 0 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
823 0 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
824 0 : rho_h = 0.0_dp
825 0 : rho_s = 0.0_dp
826 0 : drho_h = 0.0_dp
827 0 : drho_s = 0.0_dp
828 0 : 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 0 : dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
832 : END DO
833 0 : DO ir = 1, nr
834 0 : CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
835 0 : CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
836 : END DO
837 :
838 0 : IF (lsd) THEN
839 0 : IF (ispin == 1) THEN
840 0 : vxc_tmp_h(:, :, 1) = vxc_tmp_h(:, :, 1) + efac*orb_rho_set_h%rhoa(:, :, 1)
841 0 : vxc_tmp_s(:, :, 1) = vxc_tmp_s(:, :, 1) + efac*orb_rho_set_s%rhoa(:, :, 1)
842 : ELSE
843 0 : vxc_tmp_h(:, :, 2) = vxc_tmp_h(:, :, 2) + efac*orb_rho_set_h%rhob(:, :, 1)
844 0 : vxc_tmp_s(:, :, 2) = vxc_tmp_s(:, :, 2) + efac*orb_rho_set_s%rhob(:, :, 1)
845 : END IF
846 : ELSE
847 0 : vxc_tmp_h(:, :, 1) = vxc_tmp_h(:, :, 1) + efac*orb_rho_set_h%rho(:, :, 1)
848 0 : 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 0 : IF (lsd) THEN
856 0 : DO ir = 1, nr
857 0 : DO ia = 1, na
858 0 : IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) &
859 : vxc_GLLB_h(ia, ir, 1) = vxc_GLLB_h(ia, ir, 1) + &
860 0 : weight_h(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
861 0 : IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) &
862 : vxc_GLLB_h(ia, ir, 2) = vxc_GLLB_h(ia, ir, 2) + &
863 0 : weight_h(ia, ir)*vxc_tmp_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
864 0 : IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) &
865 : vxc_GLLB_s(ia, ir, 1) = vxc_GLLB_s(ia, ir, 1) + &
866 0 : weight_s(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
867 0 : IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) &
868 : vxc_GLLB_s(ia, ir, 2) = vxc_GLLB_s(ia, ir, 2) + &
869 0 : weight_s(ia, ir)*vxc_tmp_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
870 : END DO
871 : END DO
872 : ELSE
873 0 : DO ir = 1, nr
874 0 : DO ia = 1, na
875 0 : IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) &
876 : vxc_GLLB_h(ia, ir, 1) = vxc_GLLB_h(ia, ir, 1) + &
877 0 : weight_h(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
878 0 : IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) &
879 : vxc_GLLB_s(ia, ir, 1) = vxc_GLLB_s(ia, ir, 1) + &
880 0 : weight_s(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
881 : END DO
882 : END DO
883 : END IF
884 :
885 0 : vxc_SAOP_h = 0.0_dp; vxc_SAOP_s = 0.0_dp
886 :
887 0 : DO ispin = 1, nspins
888 :
889 0 : DO orb = 1, homo(ispin)
890 :
891 0 : ALLOCATE (coeff_col(nrow(ispin), 1))
892 :
893 : we_LB = EXP(-2.0_dp*(mo_eigenvalues(ispin)%array(homo(ispin)) - &
894 0 : mo_eigenvalues(ispin)%array(orb))**2)
895 0 : we_GLLB = 1.0_dp - we_LB
896 0 : IF (.NOT. lsd) THEN
897 0 : we_LB = 2.0_dp*we_LB
898 0 : we_GLLB = 2.0_dp*we_GLLB
899 : END IF
900 :
901 : vxc_tmp_h(:, :, ispin) = we_LB*vxc_LB_h(:, :, ispin) + &
902 0 : we_GLLB*vxc_GLLB_h(:, :, ispin)
903 : vxc_tmp_s(:, :, ispin) = we_LB*vxc_LB_s(:, :, ispin) + &
904 0 : we_GLLB*vxc_GLLB_s(:, :, ispin)
905 :
906 0 : CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
907 : CALL cp_fm_get_submatrix(mo_coeff(ispin)%matrix, coeff_col, &
908 0 : 1, orb, nrow(ispin), 1)
909 : CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
910 0 : 1, orb)
911 0 : CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
912 : CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
913 : matrix_v=single_mo_coeff(ispin), &
914 : ncol=ncol(ispin), &
915 0 : alpha=1.0_dp)
916 :
917 0 : DEALLOCATE (coeff_col)
918 :
919 : ! This calculates the CPC and density on the grids for every atom even though
920 : ! we need it only for iatom at the moment. It seems that to circumvent this,
921 : ! the routines must be adapted to calculate just iatom
922 : ! remap pointer
923 0 : ns = SIZE(orbital_density_matrix)
924 0 : psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
925 0 : CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
926 0 : CALL prepare_gapw_den(qs_env, local_rho_set, .FALSE.)
927 :
928 0 : rho_atom => local_rho_set%rho_atom_set(iatom)
929 0 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
930 0 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
931 0 : rho_h = 0.0_dp
932 0 : rho_s = 0.0_dp
933 0 : drho_h = 0.0_dp
934 0 : drho_s = 0.0_dp
935 0 : DO ir = 1, nr
936 : CALL calc_rho_angular(atomic_grid, harmonics, nspins, .FALSE., &
937 : ir, r_h, r_s, rho_h, rho_s, &
938 0 : dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
939 : END DO
940 0 : DO ir = 1, nr
941 0 : CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
942 0 : CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
943 : END DO
944 :
945 0 : IF (lsd) THEN
946 0 : IF (ispin == 1) THEN
947 0 : vxc_SAOP_h(:, :, 1) = vxc_SAOP_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rhoa(:, :, 1)
948 0 : vxc_SAOP_s(:, :, 1) = vxc_SAOP_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rhoa(:, :, 1)
949 : ELSE
950 0 : vxc_SAOP_h(:, :, 2) = vxc_SAOP_h(:, :, 2) + vxc_tmp_h(:, :, 2)*orb_rho_set_h%rhob(:, :, 1)
951 0 : vxc_SAOP_s(:, :, 2) = vxc_SAOP_s(:, :, 2) + vxc_tmp_s(:, :, 2)*orb_rho_set_s%rhob(:, :, 1)
952 : END IF
953 : ELSE
954 0 : vxc_SAOP_h(:, :, 1) = vxc_SAOP_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rho(:, :, 1)
955 0 : vxc_SAOP_s(:, :, 1) = vxc_SAOP_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rho(:, :, 1)
956 : END IF
957 :
958 : END DO ! orb
959 :
960 : END DO ! ispin
961 :
962 0 : IF (lsd) THEN
963 0 : DO ir = 1, nr
964 0 : DO ia = 1, na
965 0 : IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
966 0 : vxc_SAOP_h(ia, ir, 1) = vxc_SAOP_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
967 : ELSE
968 0 : vxc_SAOP_h(ia, ir, 1) = 0.0_dp
969 : END IF
970 0 : IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
971 0 : vxc_SAOP_h(ia, ir, 2) = vxc_SAOP_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
972 : ELSE
973 0 : vxc_SAOP_h(ia, ir, 2) = 0.0_dp
974 : END IF
975 0 : IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
976 0 : vxc_SAOP_s(ia, ir, 1) = vxc_SAOP_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
977 : ELSE
978 0 : vxc_SAOP_s(ia, ir, 1) = 0.0_dp
979 : END IF
980 0 : IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
981 0 : vxc_SAOP_s(ia, ir, 2) = vxc_SAOP_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
982 : ELSE
983 0 : vxc_SAOP_s(ia, ir, 2) = 0.0_dp
984 : END IF
985 : END DO
986 : END DO
987 : ELSE
988 0 : DO ir = 1, nr
989 0 : DO ia = 1, na
990 0 : IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
991 0 : vxc_SAOP_h(ia, ir, 1) = vxc_SAOP_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
992 : ELSE
993 0 : vxc_SAOP_h(ia, ir, 1) = 0.0_dp
994 : END IF
995 0 : IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
996 0 : vxc_SAOP_s(ia, ir, 1) = vxc_SAOP_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
997 : ELSE
998 0 : vxc_SAOP_s(ia, ir, 1) = 0.0_dp
999 : END IF
1000 : END DO
1001 : END DO
1002 : END IF
1003 :
1004 0 : rho_atom => rho_atom_set(iatom)
1005 0 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1006 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis, &
1007 0 : harmonics=harmonics, grid_atom=grid_atom)
1008 0 : SELECT CASE (oe_corr)
1009 : CASE (oe_lb)
1010 0 : CALL gaVxcgb_noGC(vxc_LB_h, vxc_LB_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1011 : CASE (oe_gllb)
1012 0 : CALL gaVxcgb_noGC(vxc_GLLB_h, vxc_GLLB_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1013 : CASE (oe_saop)
1014 0 : CALL gaVxcgb_noGC(vxc_SAOP_h, vxc_SAOP_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1015 : CASE default
1016 0 : CPABORT("Unknown correction!")
1017 : END SELECT
1018 :
1019 : END DO
1020 :
1021 0 : DEALLOCATE (rho_h, rho_s, weight_h, weight_s)
1022 0 : DEALLOCATE (vxc_LB_h, vxc_LB_s)
1023 0 : DEALLOCATE (vxc_GLLB_h, vxc_GLLB_s)
1024 0 : DEALLOCATE (vxc_tmp_h, vxc_tmp_s)
1025 0 : DEALLOCATE (vxc_SAOP_h, vxc_SAOP_s)
1026 0 : DEALLOCATE (drho_h, drho_s)
1027 :
1028 0 : CALL xc_dset_release(deriv_set)
1029 0 : CALL xc_rho_set_release(rho_set_h)
1030 0 : CALL xc_rho_set_release(rho_set_s)
1031 0 : CALL xc_rho_set_release(orb_rho_set_h)
1032 0 : CALL xc_rho_set_release(orb_rho_set_s)
1033 :
1034 : END DO
1035 :
1036 : ! remap pointer
1037 0 : ns = SIZE(matrix_ks)
1038 0 : ksmat(1:ns, 1:1) => matrix_ks(1:ns)
1039 0 : ns = SIZE(rho_struct_ao)
1040 0 : psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
1041 :
1042 0 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE.)
1043 :
1044 : !---------!
1045 : ! Cleanup !
1046 : !---------!
1047 0 : CALL section_vals_release(xc_fun_section_tmp)
1048 0 : CALL section_vals_release(xc_section_tmp)
1049 0 : CALL section_vals_release(xc_section_orig)
1050 :
1051 0 : CALL local_rho_set_release(local_rho_set)
1052 0 : CALL cp_fm_release(single_mo_coeff)
1053 0 : DEALLOCATE (mo_coeff, mo_eigenvalues)
1054 0 : CALL dbcsr_deallocate_matrix_set(orbital_density_matrix)
1055 :
1056 0 : END SUBROUTINE gapw_add_atomic_saop_pot
1057 :
1058 : ! **************************************************************************************************
1059 : !> \brief ...
1060 : !> \param pot ...
1061 : !> \param rho_set ...
1062 : !> \param lsd ...
1063 : !> \param spin ...
1064 : ! **************************************************************************************************
1065 20 : SUBROUTINE add_lb_pot(pot, rho_set, lsd, spin)
1066 :
1067 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pot
1068 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1069 : LOGICAL, INTENT(IN) :: lsd
1070 : INTEGER, INTENT(IN) :: spin
1071 :
1072 : REAL(KIND=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp
1073 :
1074 : INTEGER :: i, j, k
1075 : INTEGER, DIMENSION(2, 3) :: bo
1076 : REAL(KIND=dp) :: n, n_13, x, x2
1077 :
1078 200 : bo = rho_set%local_bounds
1079 :
1080 694 : DO k = bo(1, 3), bo(2, 3)
1081 24228 : DO j = bo(1, 2), bo(2, 2)
1082 447395 : DO i = bo(1, 1), bo(2, 1)
1083 446721 : IF (.NOT. lsd) THEN
1084 73875 : IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1085 73875 : n = rho_set%rho(i, j, k)/2.0_dp
1086 73875 : n_13 = n**ob3
1087 73875 : x = (rho_set%norm_drho(i, j, k)/2.0_dp)/(n*n_13)
1088 73875 : x2 = x*x
1089 73875 : pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*LOG(x + SQRT(x2 + 1.0_dp)))
1090 : END IF
1091 : ELSE
1092 349312 : IF (spin == 1) THEN
1093 174656 : IF (rho_set%rhoa(i, j, k) > rho_set%rho_cutoff) THEN
1094 174656 : n_13 = rho_set%rhoa_1_3(i, j, k)
1095 174656 : x = rho_set%norm_drhoa(i, j, k)/(rho_set%rhoa(i, j, k)*n_13)
1096 174656 : x2 = x*x
1097 174656 : pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*LOG(SQRT(x2 + 1.0_dp) + x))
1098 : END IF
1099 174656 : ELSE IF (spin == 2) THEN
1100 174656 : IF (rho_set%rhob(i, j, k) > rho_set%rho_cutoff) THEN
1101 174656 : n_13 = rho_set%rhob_1_3(i, j, k)
1102 174656 : x = rho_set%norm_drhob(i, j, k)/(rho_set%rhob(i, j, k)*n_13)
1103 174656 : x2 = x*x
1104 174656 : pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*LOG(SQRT(x2 + 1.0_dp) + x))
1105 : END IF
1106 : END IF
1107 : END IF
1108 : END DO
1109 : END DO
1110 : END DO
1111 :
1112 20 : END SUBROUTINE add_lb_pot
1113 :
1114 : ! **************************************************************************************************
1115 : !> \brief ...
1116 : !> \param pot ...
1117 : !> \param rho_set ...
1118 : !> \param e_uniform ...
1119 : !> \param lsd ...
1120 : ! **************************************************************************************************
1121 20 : SUBROUTINE calc_2excpbe(pot, rho_set, e_uniform, lsd)
1122 :
1123 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pot
1124 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1125 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: e_uniform
1126 : LOGICAL, INTENT(IN) :: lsd
1127 :
1128 : INTEGER :: i, j, k
1129 : INTEGER, DIMENSION(2, 3) :: bo
1130 : REAL(KIND=dp) :: e_unif, rho
1131 :
1132 200 : bo = rho_set%local_bounds
1133 :
1134 694 : DO k = bo(1, 3), bo(2, 3)
1135 24228 : DO j = bo(1, 2), bo(2, 2)
1136 447395 : DO i = bo(1, 1), bo(2, 1)
1137 446721 : IF (.NOT. lsd) THEN
1138 73875 : IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1139 73875 : e_unif = e_uniform(i, j, k)/rho_set%rho(i, j, k)
1140 : ELSE
1141 0 : e_unif = 0.0_dp
1142 : END IF
1143 : pot(i, j, k) = &
1144 : 2.0_dp* &
1145 : calc_ecpbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1146 : e_unif, rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1147 : 2.0_dp* &
1148 : calc_expbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1149 73875 : rho_set%rho_cutoff, rho_set%drho_cutoff)
1150 : ELSE
1151 349312 : rho = rho_set%rhoa(i, j, k) + rho_set%rhob(i, j, k)
1152 349312 : IF (rho > rho_set%rho_cutoff) THEN
1153 349312 : e_unif = e_uniform(i, j, k)/rho
1154 : ELSE
1155 0 : e_unif = 0.0_dp
1156 : END IF
1157 : pot(i, j, k) = &
1158 : 2.0_dp* &
1159 : calc_ecpbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1160 : e_unif, &
1161 : rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1162 : 2.0_dp* &
1163 : calc_expbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1164 349312 : rho_set%rho_cutoff, rho_set%drho_cutoff)
1165 : END IF
1166 : END DO
1167 : END DO
1168 : END DO
1169 :
1170 20 : END SUBROUTINE calc_2excpbe
1171 :
1172 : ! **************************************************************************************************
1173 : !> \brief ...
1174 : !> \param ra ...
1175 : !> \param rb ...
1176 : !> \param ngr ...
1177 : !> \param ec_unif ...
1178 : !> \param rc ...
1179 : !> \param ngrc ...
1180 : !> \return ...
1181 : ! **************************************************************************************************
1182 349312 : FUNCTION calc_ecpbe_u(ra, rb, ngr, ec_unif, rc, ngrc) RESULT(res)
1183 :
1184 : REAL(kind=dp), INTENT(in) :: ra, rb, ngr, ec_unif, rc, ngrc
1185 : REAL(kind=dp) :: res
1186 :
1187 : REAL(kind=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp, tb3 = 2.0_dp/3.0_dp
1188 :
1189 : REAL(kind=dp) :: A, At2, H, kf, kl, ks, phi, phi3, r, t2, &
1190 : zeta
1191 :
1192 349312 : r = ra + rb
1193 349312 : H = 0.0_dp
1194 349312 : IF (r > rc .AND. ngr > ngrc) THEN
1195 349312 : zeta = (ra - rb)/r
1196 349312 : IF (zeta > 1.0_dp) zeta = 1.0_dp ! machine precision problem
1197 349312 : IF (zeta < -1.0_dp) zeta = -1.0_dp ! machine precision problem
1198 349312 : phi = ((1.0_dp + zeta)**tb3 + (1.0_dp - zeta)**tb3)/2.0_dp
1199 349312 : phi3 = phi*phi*phi
1200 349312 : kf = (3.0_dp*r*pi*pi)**ob3
1201 349312 : ks = SQRT(4.0_dp*kf/pi)
1202 349312 : t2 = (ngr/(2.0_dp*phi*ks*r))**2
1203 349312 : A = beta_ec/gamma_saop/(EXP(-ec_unif/(gamma_saop*phi3)) - 1.0_dp)
1204 349312 : At2 = A*t2
1205 349312 : kl = (1.0_dp + At2)/(1.0_dp + At2 + At2*At2)
1206 349312 : H = gamma_saop*LOG(1.0_dp + beta_ec/gamma_saop*t2*kl)
1207 : END IF
1208 349312 : res = ec_unif + H
1209 :
1210 349312 : END FUNCTION calc_ecpbe_u
1211 :
1212 : ! **************************************************************************************************
1213 : !> \brief ...
1214 : !> \param r ...
1215 : !> \param ngr ...
1216 : !> \param ec_unif ...
1217 : !> \param rc ...
1218 : !> \param ngrc ...
1219 : !> \return ...
1220 : ! **************************************************************************************************
1221 73875 : FUNCTION calc_ecpbe_r(r, ngr, ec_unif, rc, ngrc) RESULT(res)
1222 :
1223 : REAL(kind=dp), INTENT(in) :: r, ngr, ec_unif, rc, ngrc
1224 : REAL(kind=dp) :: res
1225 :
1226 : REAL(kind=dp) :: A, At2, H, kf, kl, ks, t2
1227 :
1228 73875 : H = 0.0_dp
1229 73875 : IF (r > rc .AND. ngr > ngrc) THEN
1230 73875 : kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1231 73875 : ks = SQRT(4.0_dp*kf/pi)
1232 73875 : t2 = (ngr/(2.0_dp*ks*r))**2
1233 73875 : A = beta_ec/gamma_saop/(EXP(-ec_unif/gamma_saop) - 1.0_dp)
1234 73875 : At2 = A*t2
1235 73875 : kl = (1.0_dp + At2)/(1.0_dp + At2 + At2*At2)
1236 73875 : H = gamma_saop*LOG(1.0_dp + beta_ec/gamma_saop*t2*kl)
1237 : END IF
1238 73875 : res = ec_unif + H
1239 :
1240 73875 : END FUNCTION calc_ecpbe_r
1241 :
1242 : ! **************************************************************************************************
1243 : !> \brief ...
1244 : !> \param ra ...
1245 : !> \param rb ...
1246 : !> \param ngr ...
1247 : !> \param rc ...
1248 : !> \param ngrc ...
1249 : !> \return ...
1250 : ! **************************************************************************************************
1251 349312 : FUNCTION calc_expbe_u(ra, rb, ngr, rc, ngrc) RESULT(res)
1252 :
1253 : REAL(kind=dp), INTENT(in) :: ra, rb, ngr, rc, ngrc
1254 : REAL(kind=dp) :: res
1255 :
1256 : REAL(kind=dp) :: r
1257 :
1258 349312 : r = ra + rb
1259 349312 : res = calc_expbe_r(r, ngr, rc, ngrc)
1260 :
1261 349312 : END FUNCTION calc_expbe_u
1262 :
1263 : ! **************************************************************************************************
1264 : !> \brief ...
1265 : !> \param r ...
1266 : !> \param ngr ...
1267 : !> \param rc ...
1268 : !> \param ngrc ...
1269 : !> \return ...
1270 : ! **************************************************************************************************
1271 423187 : FUNCTION calc_expbe_r(r, ngr, rc, ngrc) RESULT(res)
1272 :
1273 : REAL(kind=dp), INTENT(in) :: r, ngr, rc, ngrc
1274 : REAL(kind=dp) :: res
1275 :
1276 : REAL(kind=dp) :: ex_unif, fx, kf, s
1277 :
1278 423187 : IF (r > rc) THEN
1279 423187 : kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1280 423187 : ex_unif = -3.0_dp*kf/(4.0_dp*pi)
1281 423187 : fx = 1.0_dp
1282 423187 : IF (ngr > ngrc) THEN
1283 423187 : s = ngr/(2.0_dp*kf*r)
1284 423187 : fx = fx + kappa - kappa/(1.0_dp + mu*s*s/kappa)
1285 : END IF
1286 423187 : res = ex_unif*fx
1287 : ELSE
1288 : res = 0.0_dp
1289 : END IF
1290 :
1291 423187 : END FUNCTION calc_expbe_r
1292 :
1293 : END MODULE xc_pot_saop
|
