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