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