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 : 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 : 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) THEN
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 : END IF
862 0 : IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
863 : vxc_GLLB_h(ia, ir, 2) = vxc_GLLB_h(ia, ir, 2) + &
864 0 : weight_h(ia, ir)*vxc_tmp_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
865 : END IF
866 0 : IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
867 : vxc_GLLB_s(ia, ir, 1) = vxc_GLLB_s(ia, ir, 1) + &
868 0 : weight_s(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
869 : END IF
870 0 : IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
871 : vxc_GLLB_s(ia, ir, 2) = vxc_GLLB_s(ia, ir, 2) + &
872 0 : weight_s(ia, ir)*vxc_tmp_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
873 : END IF
874 : END DO
875 : END DO
876 : ELSE
877 0 : DO ir = 1, nr
878 0 : DO ia = 1, na
879 0 : IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
880 : vxc_GLLB_h(ia, ir, 1) = vxc_GLLB_h(ia, ir, 1) + &
881 0 : weight_h(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
882 : END IF
883 0 : IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
884 : vxc_GLLB_s(ia, ir, 1) = vxc_GLLB_s(ia, ir, 1) + &
885 0 : weight_s(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
886 : END IF
887 : END DO
888 : END DO
889 : END IF
890 :
891 0 : vxc_SAOP_h = 0.0_dp; vxc_SAOP_s = 0.0_dp
892 :
893 0 : DO ispin = 1, nspins
894 :
895 0 : DO orb = 1, homo(ispin)
896 :
897 0 : ALLOCATE (coeff_col(nrow(ispin), 1))
898 :
899 : we_LB = EXP(-2.0_dp*(mo_eigenvalues(ispin)%array(homo(ispin)) - &
900 0 : mo_eigenvalues(ispin)%array(orb))**2)
901 0 : we_GLLB = 1.0_dp - we_LB
902 0 : IF (.NOT. lsd) THEN
903 0 : we_LB = 2.0_dp*we_LB
904 0 : we_GLLB = 2.0_dp*we_GLLB
905 : END IF
906 :
907 : vxc_tmp_h(:, :, ispin) = we_LB*vxc_LB_h(:, :, ispin) + &
908 0 : we_GLLB*vxc_GLLB_h(:, :, ispin)
909 : vxc_tmp_s(:, :, ispin) = we_LB*vxc_LB_s(:, :, ispin) + &
910 0 : we_GLLB*vxc_GLLB_s(:, :, ispin)
911 :
912 0 : CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
913 : CALL cp_fm_get_submatrix(mo_coeff(ispin)%matrix, coeff_col, &
914 0 : 1, orb, nrow(ispin), 1)
915 : CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
916 0 : 1, orb)
917 0 : CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
918 : CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
919 : matrix_v=single_mo_coeff(ispin), &
920 : ncol=ncol(ispin), &
921 0 : alpha=1.0_dp)
922 :
923 0 : DEALLOCATE (coeff_col)
924 :
925 : ! This calculates the CPC and density on the grids for every atom even though
926 : ! we need it only for iatom at the moment. It seems that to circumvent this,
927 : ! the routines must be adapted to calculate just iatom
928 : ! remap pointer
929 0 : ns = SIZE(orbital_density_matrix)
930 0 : psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
931 0 : CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
932 0 : CALL prepare_gapw_den(qs_env, local_rho_set, .FALSE.)
933 :
934 0 : rho_atom => local_rho_set%rho_atom_set(iatom)
935 0 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
936 0 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
937 0 : rho_h = 0.0_dp
938 0 : rho_s = 0.0_dp
939 0 : drho_h = 0.0_dp
940 0 : drho_s = 0.0_dp
941 0 : DO ir = 1, nr
942 : CALL calc_rho_angular(atomic_grid, harmonics, nspins, .FALSE., &
943 : ir, r_h, r_s, rho_h, rho_s, &
944 0 : dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
945 : END DO
946 0 : DO ir = 1, nr
947 0 : CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
948 0 : CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
949 : END DO
950 :
951 0 : IF (lsd) THEN
952 0 : IF (ispin == 1) THEN
953 0 : vxc_SAOP_h(:, :, 1) = vxc_SAOP_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rhoa(:, :, 1)
954 0 : vxc_SAOP_s(:, :, 1) = vxc_SAOP_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rhoa(:, :, 1)
955 : ELSE
956 0 : vxc_SAOP_h(:, :, 2) = vxc_SAOP_h(:, :, 2) + vxc_tmp_h(:, :, 2)*orb_rho_set_h%rhob(:, :, 1)
957 0 : vxc_SAOP_s(:, :, 2) = vxc_SAOP_s(:, :, 2) + vxc_tmp_s(:, :, 2)*orb_rho_set_s%rhob(:, :, 1)
958 : END IF
959 : ELSE
960 0 : vxc_SAOP_h(:, :, 1) = vxc_SAOP_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rho(:, :, 1)
961 0 : vxc_SAOP_s(:, :, 1) = vxc_SAOP_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rho(:, :, 1)
962 : END IF
963 :
964 : END DO ! orb
965 :
966 : END DO ! ispin
967 :
968 0 : IF (lsd) THEN
969 0 : DO ir = 1, nr
970 0 : DO ia = 1, na
971 0 : IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
972 0 : vxc_SAOP_h(ia, ir, 1) = vxc_SAOP_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
973 : ELSE
974 0 : vxc_SAOP_h(ia, ir, 1) = 0.0_dp
975 : END IF
976 0 : IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
977 0 : vxc_SAOP_h(ia, ir, 2) = vxc_SAOP_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
978 : ELSE
979 0 : vxc_SAOP_h(ia, ir, 2) = 0.0_dp
980 : END IF
981 0 : IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
982 0 : vxc_SAOP_s(ia, ir, 1) = vxc_SAOP_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
983 : ELSE
984 0 : vxc_SAOP_s(ia, ir, 1) = 0.0_dp
985 : END IF
986 0 : IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
987 0 : vxc_SAOP_s(ia, ir, 2) = vxc_SAOP_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
988 : ELSE
989 0 : vxc_SAOP_s(ia, ir, 2) = 0.0_dp
990 : END IF
991 : END DO
992 : END DO
993 : ELSE
994 0 : DO ir = 1, nr
995 0 : DO ia = 1, na
996 0 : IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
997 0 : vxc_SAOP_h(ia, ir, 1) = vxc_SAOP_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
998 : ELSE
999 0 : vxc_SAOP_h(ia, ir, 1) = 0.0_dp
1000 : END IF
1001 0 : IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
1002 0 : vxc_SAOP_s(ia, ir, 1) = vxc_SAOP_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
1003 : ELSE
1004 0 : vxc_SAOP_s(ia, ir, 1) = 0.0_dp
1005 : END IF
1006 : END DO
1007 : END DO
1008 : END IF
1009 :
1010 0 : rho_atom => rho_atom_set(iatom)
1011 0 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1012 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis, &
1013 0 : harmonics=harmonics, grid_atom=grid_atom)
1014 0 : SELECT CASE (oe_corr)
1015 : CASE (oe_lb)
1016 0 : CALL gaVxcgb_noGC(vxc_LB_h, vxc_LB_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1017 : CASE (oe_gllb)
1018 0 : CALL gaVxcgb_noGC(vxc_GLLB_h, vxc_GLLB_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1019 : CASE (oe_saop)
1020 0 : CALL gaVxcgb_noGC(vxc_SAOP_h, vxc_SAOP_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1021 : CASE default
1022 0 : CPABORT("Unknown correction!")
1023 : END SELECT
1024 :
1025 : END DO
1026 :
1027 0 : DEALLOCATE (rho_h, rho_s, weight_h, weight_s)
1028 0 : DEALLOCATE (vxc_LB_h, vxc_LB_s)
1029 0 : DEALLOCATE (vxc_GLLB_h, vxc_GLLB_s)
1030 0 : DEALLOCATE (vxc_tmp_h, vxc_tmp_s)
1031 0 : DEALLOCATE (vxc_SAOP_h, vxc_SAOP_s)
1032 0 : DEALLOCATE (drho_h, drho_s)
1033 :
1034 0 : CALL xc_dset_release(deriv_set)
1035 0 : CALL xc_rho_set_release(rho_set_h)
1036 0 : CALL xc_rho_set_release(rho_set_s)
1037 0 : CALL xc_rho_set_release(orb_rho_set_h)
1038 0 : CALL xc_rho_set_release(orb_rho_set_s)
1039 :
1040 : END DO
1041 :
1042 : ! remap pointer
1043 0 : ns = SIZE(matrix_ks)
1044 0 : ksmat(1:ns, 1:1) => matrix_ks(1:ns)
1045 0 : ns = SIZE(rho_struct_ao)
1046 0 : psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
1047 :
1048 0 : CALL update_ks_atom(qs_env, ksmat, psmat, forces=.FALSE.)
1049 :
1050 : !---------!
1051 : ! Cleanup !
1052 : !---------!
1053 0 : CALL section_vals_release(xc_fun_section_tmp)
1054 0 : CALL section_vals_release(xc_section_tmp)
1055 0 : CALL section_vals_release(xc_section_orig)
1056 :
1057 0 : CALL local_rho_set_release(local_rho_set)
1058 0 : CALL cp_fm_release(single_mo_coeff)
1059 0 : DEALLOCATE (mo_coeff, mo_eigenvalues)
1060 0 : CALL dbcsr_deallocate_matrix_set(orbital_density_matrix)
1061 :
1062 0 : END SUBROUTINE gapw_add_atomic_saop_pot
1063 :
1064 : ! **************************************************************************************************
1065 : !> \brief ...
1066 : !> \param pot ...
1067 : !> \param rho_set ...
1068 : !> \param lsd ...
1069 : !> \param spin ...
1070 : ! **************************************************************************************************
1071 20 : SUBROUTINE add_lb_pot(pot, rho_set, lsd, spin)
1072 :
1073 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pot
1074 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1075 : LOGICAL, INTENT(IN) :: lsd
1076 : INTEGER, INTENT(IN) :: spin
1077 :
1078 : REAL(KIND=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp
1079 :
1080 : INTEGER :: i, j, k
1081 : INTEGER, DIMENSION(2, 3) :: bo
1082 : REAL(KIND=dp) :: n, n_13, x, x2
1083 :
1084 200 : bo = rho_set%local_bounds
1085 :
1086 694 : DO k = bo(1, 3), bo(2, 3)
1087 24228 : DO j = bo(1, 2), bo(2, 2)
1088 447395 : DO i = bo(1, 1), bo(2, 1)
1089 446721 : IF (.NOT. lsd) THEN
1090 73875 : IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1091 73875 : n = rho_set%rho(i, j, k)/2.0_dp
1092 73875 : n_13 = n**ob3
1093 73875 : x = (rho_set%norm_drho(i, j, k)/2.0_dp)/(n*n_13)
1094 73875 : x2 = x*x
1095 73875 : pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*LOG(x + SQRT(x2 + 1.0_dp)))
1096 : END IF
1097 : ELSE
1098 349312 : IF (spin == 1) THEN
1099 174656 : IF (rho_set%rhoa(i, j, k) > rho_set%rho_cutoff) THEN
1100 174656 : n_13 = rho_set%rhoa_1_3(i, j, k)
1101 174656 : x = rho_set%norm_drhoa(i, j, k)/(rho_set%rhoa(i, j, k)*n_13)
1102 174656 : x2 = x*x
1103 174656 : pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*LOG(SQRT(x2 + 1.0_dp) + x))
1104 : END IF
1105 174656 : ELSE IF (spin == 2) THEN
1106 174656 : IF (rho_set%rhob(i, j, k) > rho_set%rho_cutoff) THEN
1107 174656 : n_13 = rho_set%rhob_1_3(i, j, k)
1108 174656 : x = rho_set%norm_drhob(i, j, k)/(rho_set%rhob(i, j, k)*n_13)
1109 174656 : x2 = x*x
1110 174656 : pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*LOG(SQRT(x2 + 1.0_dp) + x))
1111 : END IF
1112 : END IF
1113 : END IF
1114 : END DO
1115 : END DO
1116 : END DO
1117 :
1118 20 : END SUBROUTINE add_lb_pot
1119 :
1120 : ! **************************************************************************************************
1121 : !> \brief ...
1122 : !> \param pot ...
1123 : !> \param rho_set ...
1124 : !> \param e_uniform ...
1125 : !> \param lsd ...
1126 : ! **************************************************************************************************
1127 20 : SUBROUTINE calc_2excpbe(pot, rho_set, e_uniform, lsd)
1128 :
1129 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: pot
1130 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1131 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: e_uniform
1132 : LOGICAL, INTENT(IN) :: lsd
1133 :
1134 : INTEGER :: i, j, k
1135 : INTEGER, DIMENSION(2, 3) :: bo
1136 : REAL(KIND=dp) :: e_unif, rho
1137 :
1138 200 : bo = rho_set%local_bounds
1139 :
1140 694 : DO k = bo(1, 3), bo(2, 3)
1141 24228 : DO j = bo(1, 2), bo(2, 2)
1142 447395 : DO i = bo(1, 1), bo(2, 1)
1143 446721 : IF (.NOT. lsd) THEN
1144 73875 : IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1145 73875 : e_unif = e_uniform(i, j, k)/rho_set%rho(i, j, k)
1146 : ELSE
1147 0 : e_unif = 0.0_dp
1148 : END IF
1149 : pot(i, j, k) = &
1150 : 2.0_dp* &
1151 : calc_ecpbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1152 : e_unif, rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1153 : 2.0_dp* &
1154 : calc_expbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1155 73875 : rho_set%rho_cutoff, rho_set%drho_cutoff)
1156 : ELSE
1157 349312 : rho = rho_set%rhoa(i, j, k) + rho_set%rhob(i, j, k)
1158 349312 : IF (rho > rho_set%rho_cutoff) THEN
1159 349312 : e_unif = e_uniform(i, j, k)/rho
1160 : ELSE
1161 0 : e_unif = 0.0_dp
1162 : END IF
1163 : pot(i, j, k) = &
1164 : 2.0_dp* &
1165 : calc_ecpbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1166 : e_unif, &
1167 : rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1168 : 2.0_dp* &
1169 : calc_expbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1170 349312 : rho_set%rho_cutoff, rho_set%drho_cutoff)
1171 : END IF
1172 : END DO
1173 : END DO
1174 : END DO
1175 :
1176 20 : END SUBROUTINE calc_2excpbe
1177 :
1178 : ! **************************************************************************************************
1179 : !> \brief ...
1180 : !> \param ra ...
1181 : !> \param rb ...
1182 : !> \param ngr ...
1183 : !> \param ec_unif ...
1184 : !> \param rc ...
1185 : !> \param ngrc ...
1186 : !> \return ...
1187 : ! **************************************************************************************************
1188 349312 : FUNCTION calc_ecpbe_u(ra, rb, ngr, ec_unif, rc, ngrc) RESULT(res)
1189 :
1190 : REAL(kind=dp), INTENT(in) :: ra, rb, ngr, ec_unif, rc, ngrc
1191 : REAL(kind=dp) :: res
1192 :
1193 : REAL(kind=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp, tb3 = 2.0_dp/3.0_dp
1194 :
1195 : REAL(kind=dp) :: A, At2, H, kf, kl, ks, phi, phi3, r, t2, &
1196 : zeta
1197 :
1198 349312 : r = ra + rb
1199 349312 : H = 0.0_dp
1200 349312 : IF (r > rc .AND. ngr > ngrc) THEN
1201 349312 : zeta = (ra - rb)/r
1202 349312 : IF (zeta > 1.0_dp) zeta = 1.0_dp ! machine precision problem
1203 : IF (zeta < -1.0_dp) zeta = -1.0_dp ! machine precision problem
1204 349312 : phi = ((1.0_dp + zeta)**tb3 + (1.0_dp - zeta)**tb3)/2.0_dp
1205 349312 : phi3 = phi*phi*phi
1206 349312 : kf = (3.0_dp*r*pi*pi)**ob3
1207 349312 : ks = SQRT(4.0_dp*kf/pi)
1208 349312 : t2 = (ngr/(2.0_dp*phi*ks*r))**2
1209 349312 : A = beta_ec/gamma_saop/(EXP(-ec_unif/(gamma_saop*phi3)) - 1.0_dp)
1210 349312 : At2 = A*t2
1211 349312 : kl = (1.0_dp + At2)/(1.0_dp + At2 + At2*At2)
1212 349312 : H = gamma_saop*LOG(1.0_dp + beta_ec/gamma_saop*t2*kl)
1213 : END IF
1214 349312 : res = ec_unif + H
1215 :
1216 349312 : END FUNCTION calc_ecpbe_u
1217 :
1218 : ! **************************************************************************************************
1219 : !> \brief ...
1220 : !> \param r ...
1221 : !> \param ngr ...
1222 : !> \param ec_unif ...
1223 : !> \param rc ...
1224 : !> \param ngrc ...
1225 : !> \return ...
1226 : ! **************************************************************************************************
1227 73875 : FUNCTION calc_ecpbe_r(r, ngr, ec_unif, rc, ngrc) RESULT(res)
1228 :
1229 : REAL(kind=dp), INTENT(in) :: r, ngr, ec_unif, rc, ngrc
1230 : REAL(kind=dp) :: res
1231 :
1232 : REAL(kind=dp) :: A, At2, H, kf, kl, ks, t2
1233 :
1234 73875 : H = 0.0_dp
1235 73875 : IF (r > rc .AND. ngr > ngrc) THEN
1236 73875 : kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1237 73875 : ks = SQRT(4.0_dp*kf/pi)
1238 73875 : t2 = (ngr/(2.0_dp*ks*r))**2
1239 73875 : A = beta_ec/gamma_saop/(EXP(-ec_unif/gamma_saop) - 1.0_dp)
1240 73875 : At2 = A*t2
1241 73875 : kl = (1.0_dp + At2)/(1.0_dp + At2 + At2*At2)
1242 73875 : H = gamma_saop*LOG(1.0_dp + beta_ec/gamma_saop*t2*kl)
1243 : END IF
1244 73875 : res = ec_unif + H
1245 :
1246 73875 : END FUNCTION calc_ecpbe_r
1247 :
1248 : ! **************************************************************************************************
1249 : !> \brief ...
1250 : !> \param ra ...
1251 : !> \param rb ...
1252 : !> \param ngr ...
1253 : !> \param rc ...
1254 : !> \param ngrc ...
1255 : !> \return ...
1256 : ! **************************************************************************************************
1257 349312 : FUNCTION calc_expbe_u(ra, rb, ngr, rc, ngrc) RESULT(res)
1258 :
1259 : REAL(kind=dp), INTENT(in) :: ra, rb, ngr, rc, ngrc
1260 : REAL(kind=dp) :: res
1261 :
1262 : REAL(kind=dp) :: r
1263 :
1264 349312 : r = ra + rb
1265 349312 : res = calc_expbe_r(r, ngr, rc, ngrc)
1266 :
1267 349312 : END FUNCTION calc_expbe_u
1268 :
1269 : ! **************************************************************************************************
1270 : !> \brief ...
1271 : !> \param r ...
1272 : !> \param ngr ...
1273 : !> \param rc ...
1274 : !> \param ngrc ...
1275 : !> \return ...
1276 : ! **************************************************************************************************
1277 423187 : FUNCTION calc_expbe_r(r, ngr, rc, ngrc) RESULT(res)
1278 :
1279 : REAL(kind=dp), INTENT(in) :: r, ngr, rc, ngrc
1280 : REAL(kind=dp) :: res
1281 :
1282 : REAL(kind=dp) :: ex_unif, fx, kf, s
1283 :
1284 423187 : IF (r > rc) THEN
1285 423187 : kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1286 423187 : ex_unif = -3.0_dp*kf/(4.0_dp*pi)
1287 423187 : fx = 1.0_dp
1288 423187 : IF (ngr > ngrc) THEN
1289 423187 : s = ngr/(2.0_dp*kf*r)
1290 423187 : fx = fx + kappa - kappa/(1.0_dp + mu*s*s/kappa)
1291 : END IF
1292 423187 : res = ex_unif*fx
1293 : ELSE
1294 : res = 0.0_dp
1295 : END IF
1296 :
1297 423187 : END FUNCTION calc_expbe_r
1298 :
1299 : END MODULE xc_pot_saop
|
