Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief This functional is a combination of three different exchange hole
10 : !> models. The ingredients are:
11 : !>
12 : !> 1. Becke Roussel exchange hole
13 : !> 2. PBE exchange hole
14 : !> 3. LDA exchange hole
15 : !>
16 : !> The full functionals is given as follows
17 : !>
18 : !> Fx = eps_lr_lda/eps_lr_br
19 : !> Fcorr = alpha/( exp( (Fx-mu)/N ) + 1)
20 : !> rhox = Fcorr * eps_lr_pbe + (1-Fcorr)*eps_lr_br
21 : !> eps = int_{R}^{\infty} rhox*s*ds
22 : !>
23 : !> with alpha, mu and N fitting parameters
24 : !> \par History
25 : !> 01.2009 created [mguidon]
26 : !> \author mguidon
27 : ! **************************************************************************************************
28 :
29 : MODULE xc_xbr_pbe_lda_hole_t_c_lr
30 :
31 : USE input_section_types, ONLY: section_vals_type,&
32 : section_vals_val_get
33 : USE kinds, ONLY: dp
34 : USE mathconstants, ONLY: pi
35 : USE xc_derivative_desc, ONLY: &
36 : deriv_laplace_rho, deriv_laplace_rhoa, deriv_laplace_rhob, deriv_norm_drho, &
37 : deriv_norm_drhoa, deriv_norm_drhob, deriv_rho, deriv_rhoa, deriv_rhob, deriv_tau, &
38 : deriv_tau_a, deriv_tau_b
39 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
40 : xc_dset_get_derivative
41 : USE xc_derivative_types, ONLY: xc_derivative_get,&
42 : xc_derivative_type
43 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
44 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
45 : xc_rho_set_type
46 : USE xc_xbecke_roussel, ONLY: x_br_lsd_y_gt_0,&
47 : x_br_lsd_y_gt_0_cutoff,&
48 : x_br_lsd_y_lte_0,&
49 : x_br_lsd_y_lte_0_cutoff
50 : USE xc_xlda_hole_t_c_lr, ONLY: xlda_hole_t_c_lr_lda_calc_0
51 : USE xc_xpbe_hole_t_c_lr, ONLY: xpbe_hole_t_c_lr_lda_calc_1,&
52 : xpbe_hole_t_c_lr_lda_calc_2
53 : #include "../base/base_uses.f90"
54 :
55 : IMPLICIT NONE
56 : PRIVATE
57 :
58 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbr_pbe_lda_hole_t_c_lr'
60 :
61 : REAL(dp), PARAMETER, PRIVATE :: br_a1 = 1.5255251812009530_dp, &
62 : br_a2 = 0.4576575543602858_dp, &
63 : br_a3 = 0.4292036732051034_dp, &
64 : br_c0 = 0.7566445420735584_dp, &
65 : br_c1 = -2.6363977871370960_dp, &
66 : br_c2 = 5.4745159964232880_dp, &
67 : br_c3 = -12.657308127108290_dp, &
68 : br_c4 = 4.1250584725121360_dp, &
69 : br_c5 = -30.425133957163840_dp, &
70 : br_b0 = 0.4771976183772063_dp, &
71 : br_b1 = -1.7799813494556270_dp, &
72 : br_b2 = 3.8433841862302150_dp, &
73 : br_b3 = -9.5912050880518490_dp, &
74 : br_b4 = 2.1730180285916720_dp, &
75 : br_b5 = -30.425133851603660_dp, &
76 : br_d0 = 0.00004435009886795587_dp, &
77 : br_d1 = 0.58128653604457910_dp, &
78 : br_d2 = 66.742764515940610_dp, &
79 : br_d3 = 434.26780897229770_dp, &
80 : br_d4 = 824.7765766052239000_dp, &
81 : br_d5 = 1657.9652731582120_dp, &
82 : br_e0 = 0.00003347285060926091_dp, &
83 : br_e1 = 0.47917931023971350_dp, &
84 : br_e2 = 62.392268338574240_dp, &
85 : br_e3 = 463.14816427938120_dp, &
86 : br_e4 = 785.2360350104029000_dp, &
87 : br_e5 = 1657.962968223273000000_dp, &
88 : br_BB = 2.085749716493756_dp
89 :
90 : REAL(dp), PARAMETER, PRIVATE :: smax = 8.572844_dp, &
91 : scutoff = 8.3_dp, &
92 : sconst = 18.79622316_dp, &
93 : gcutoff = 0.08_dp
94 :
95 : REAL(dp), PARAMETER, PRIVATE :: alpha = 0.3956891_dp, &
96 : N = -0.0009800242_dp, &
97 : mu = 0.00118684_dp
98 :
99 : PUBLIC :: xbr_pbe_lda_hole_tc_lr_lda_info, &
100 : xbr_pbe_lda_hole_tc_lr_lsd_info, &
101 : xbr_pbe_lda_hole_tc_lr_lda_eval, &
102 : xbr_pbe_lda_hole_tc_lr_lsd_eval
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief return various information on the functional
107 : !> \param reference string with the reference of the actual functional
108 : !> \param shortform string with the shortform of the functional name
109 : !> \param needs the components needed by this functional are set to
110 : !> true (does not set the unneeded components to false)
111 : !> \param max_deriv ...
112 : !> \author mguidon (01.2009)
113 : ! **************************************************************************************************
114 205 : SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_info(reference, shortform, needs, max_deriv)
115 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
116 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
117 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
118 :
119 205 : IF (PRESENT(reference)) THEN
120 1 : reference = "{LDA version}"
121 : END IF
122 205 : IF (PRESENT(shortform)) THEN
123 1 : shortform = "{LDA}"
124 : END IF
125 :
126 205 : IF (PRESENT(needs)) THEN
127 204 : needs%rho = .TRUE.
128 204 : needs%norm_drho = .TRUE.
129 204 : needs%tau = .TRUE.
130 204 : needs%laplace_rho = .TRUE.
131 : END IF
132 :
133 205 : IF (PRESENT(max_deriv)) max_deriv = 1
134 :
135 205 : END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_info
136 :
137 : ! **************************************************************************************************
138 : !> \brief return various information on the functional
139 : !> \param reference string with the reference of the actual functional
140 : !> \param shortform string with the shortform of the functional name
141 : !> \param needs the components needed by this functional are set to
142 : !> true (does not set the unneeded components to false)
143 : !> \param max_deriv ...
144 : !> \author mguidon (01.2009)
145 : ! **************************************************************************************************
146 63 : SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_info(reference, shortform, needs, max_deriv)
147 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
148 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
149 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
150 :
151 63 : IF (PRESENT(reference)) THEN
152 1 : reference = "{LDA version}"
153 : END IF
154 63 : IF (PRESENT(shortform)) THEN
155 1 : shortform = "{LDA}"
156 : END IF
157 :
158 63 : IF (PRESENT(needs)) THEN
159 62 : needs%rho_spin = .TRUE.
160 62 : needs%norm_drho_spin = .TRUE.
161 62 : needs%tau_spin = .TRUE.
162 62 : needs%laplace_rho_spin = .TRUE.
163 : END IF
164 63 : IF (PRESENT(max_deriv)) max_deriv = 1
165 :
166 63 : END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_info
167 :
168 : ! **************************************************************************************************
169 : !> \brief Intermediate routine that gets grids, derivatives and some params
170 : !> \param rho_set the density where you want to evaluate the functional
171 : !> \param deriv_set place where to store the functional derivatives (they are
172 : !> added to the derivatives)
173 : !> \param grad_deriv degree of the derivative that should be evaluated,
174 : !> if positive all the derivatives up to the given degree are evaluated,
175 : !> if negative only the given degree is calculated
176 : !> \param params parameters for functional
177 : !> \author mguidon (01.2009)
178 : ! **************************************************************************************************
179 1000 : SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_eval(rho_set, deriv_set, grad_deriv, params)
180 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
181 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
182 : INTEGER, INTENT(in) :: grad_deriv
183 : TYPE(section_vals_type), POINTER :: params
184 :
185 : CHARACTER(len=*), PARAMETER :: routineN = 'xbr_pbe_lda_hole_tc_lr_lda_eval'
186 :
187 : INTEGER :: handle, npoints
188 : INTEGER, DIMENSION(2, 3) :: bo
189 : REAL(dp) :: gamma, R, sx
190 : REAL(kind=dp) :: epsilon_rho
191 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
192 200 : POINTER :: dummy, e_0, e_laplace_rho, e_ndrho, &
193 200 : e_rho, e_tau, laplace_rho, norm_drho, &
194 200 : rho, tau
195 : TYPE(xc_derivative_type), POINTER :: deriv
196 :
197 200 : CALL timeset(routineN, handle)
198 :
199 : CALL xc_rho_set_get(rho_set, rho=rho, norm_drho=norm_drho, &
200 : tau=tau, laplace_rho=laplace_rho, local_bounds=bo, &
201 200 : rho_cutoff=epsilon_rho)
202 200 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
203 :
204 200 : dummy => rho
205 :
206 200 : e_0 => dummy
207 200 : e_rho => dummy
208 200 : e_ndrho => dummy
209 200 : e_tau => dummy
210 200 : e_laplace_rho => dummy
211 :
212 200 : IF (grad_deriv >= 0) THEN
213 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
214 200 : allocate_deriv=.TRUE.)
215 200 : CALL xc_derivative_get(deriv, deriv_data=e_0)
216 : END IF
217 200 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
218 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
219 102 : allocate_deriv=.TRUE.)
220 102 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
221 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
222 102 : allocate_deriv=.TRUE.)
223 102 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
224 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
225 102 : allocate_deriv=.TRUE.)
226 102 : CALL xc_derivative_get(deriv, deriv_data=e_tau)
227 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rho], &
228 102 : allocate_deriv=.TRUE.)
229 102 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rho)
230 : END IF
231 200 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
232 0 : CPABORT("derivatives bigger than 1 not implemented")
233 : END IF
234 :
235 200 : CALL section_vals_val_get(params, "scale_x", r_val=sx)
236 200 : CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
237 200 : CALL section_vals_val_get(params, "GAMMA", r_val=gamma)
238 :
239 200 : IF (R == 0.0_dp) THEN
240 0 : CPABORT("Cutoff_Radius 0.0 not implemented")
241 : END IF
242 :
243 : !$OMP PARALLEL DEFAULT(NONE) &
244 : !$OMP SHARED(rho, norm_drho, laplace_rho, tau, e_0, e_rho) &
245 : !$OMP SHARED(e_ndrho, e_tau, e_laplace_rho, grad_deriv) &
246 : !$OMP SHARED(npoints, epsilon_rho) &
247 200 : !$OMP SHARED(sx, r, gamma)
248 :
249 : CALL xbr_pbe_lda_hole_tc_lr_lda_calc(rho=rho, norm_drho=norm_drho, &
250 : laplace_rho=laplace_rho, tau=tau, e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
251 : e_tau=e_tau, e_laplace_rho=e_laplace_rho, grad_deriv=grad_deriv, &
252 : npoints=npoints, epsilon_rho=epsilon_rho, sx=sx, R=R, gamma=gamma)
253 :
254 : !$OMP END PARALLEL
255 :
256 200 : CALL timestop(handle)
257 200 : END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_eval
258 :
259 : ! **************************************************************************************************
260 : !> \brief Low level routine that calls the three involved holes and puts them
261 : !> together
262 : !> \param rho values on the grid
263 : !> \param norm_drho values on the grid
264 : !> \param laplace_rho values on the grid
265 : !> \param tau values on the grid
266 : !> \param e_0 derivatives on the grid
267 : !> \param e_rho derivatives on the grid
268 : !> \param e_ndrho derivatives on the grid
269 : !> \param e_tau derivatives on the grid
270 : !> \param e_laplace_rho derivatives on the grid
271 : !> \param grad_deriv degree of the derivative that should be evaluated,
272 : !> if positive all the derivatives up to the given degree are evaluated,
273 : !> if negative only the given degree is calculated
274 : !> \param npoints number of gridpoints
275 : !> \param epsilon_rho cutoffs
276 : !> \param sx parameters for functional
277 : !> \param R parameters for functional
278 : !> \param gamma parameters for functional
279 : !> \author mguidon (01.2009)
280 : ! **************************************************************************************************
281 200 : SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
282 200 : e_ndrho, e_tau, e_laplace_rho, grad_deriv, npoints, &
283 : epsilon_rho, sx, R, gamma)
284 :
285 : INTEGER, INTENT(in) :: npoints, grad_deriv
286 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_laplace_rho, e_tau, e_ndrho, e_rho, e_0
287 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: tau, laplace_rho, norm_drho, rho
288 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R, gamma
289 :
290 : INTEGER :: ip
291 : REAL(dp) :: dFermi_dlaplace_rho, dFermi_dndrho, dFermi_drho, dFermi_dtau, e_0_br, e_0_lda, &
292 : e_0_pbe, e_laplace_rho_br, e_ndrho_br, e_ndrho_pbe, e_rho_br, e_rho_lda, e_rho_pbe, &
293 : e_tau_br, Fermi, Fx, my_laplace_rho, my_ndrho, my_rho, my_tau, ss, ss2, sscale, t1, t15, &
294 : t16, t2, t3, t4, t5, t6, t7, t8, t9, yval
295 :
296 : !$OMP DO
297 :
298 : DO ip = 1, npoints
299 6400000 : my_rho = 0.5_dp*MAX(rho(ip), 0.0_dp)
300 6400000 : IF (my_rho > epsilon_rho) THEN
301 6400000 : my_ndrho = 0.5_dp*MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp)
302 6400000 : my_tau = 0.5_dp*MAX(EPSILON(0.0_dp)*1.e4_dp, tau(ip))
303 6400000 : my_laplace_rho = 0.5_dp*laplace_rho(ip)
304 :
305 : ! ** We calculate first the Becke-Roussel part, saving everything in local variables
306 6400000 : t1 = pi**(0.1e1_dp/0.3e1_dp)
307 6400000 : t2 = t1**2
308 6400000 : t3 = my_rho**(0.1e1_dp/0.3e1_dp)
309 6400000 : t4 = t3**2
310 6400000 : t5 = t4*my_rho
311 6400000 : t8 = my_ndrho**2
312 6400000 : t9 = 0.1e1_dp/my_rho
313 6400000 : t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
314 6400000 : t16 = 0.1e1_dp/t15
315 6400000 : yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
316 :
317 6400000 : e_0_br = 0.0_dp
318 6400000 : e_rho_br = 0.0_dp
319 6400000 : e_ndrho_br = 0.0_dp
320 6400000 : e_tau_br = 0.0_dp
321 6400000 : e_laplace_rho_br = 0.0_dp
322 :
323 6400000 : IF (R == 0.0_dp) THEN
324 0 : IF (yval <= 0.0_dp) THEN
325 : CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
326 : e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
327 0 : sx, gamma, grad_deriv)
328 : ! VERY UGLY HACK e_0 has to multiplied by the factor 2
329 0 : e_0_br = 2.0_dp*e_0_br
330 : ELSE
331 : CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
332 : e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
333 0 : sx, gamma, grad_deriv)
334 : ! VERY UGLY HACK e_0 has to multiplied by the factor 2
335 0 : e_0_br = 2.0_dp*e_0_br
336 : END IF
337 : ELSE
338 6400000 : IF (yval <= 0.0_dp) THEN
339 : CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
340 : e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
341 143257 : sx, R, gamma, grad_deriv)
342 : ! VERY UGLY HACK e_0 has to multiplied by the factor 2
343 143257 : e_0_br = 2.0_dp*e_0_br
344 : ELSE
345 : CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
346 : e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
347 6256743 : sx, R, gamma, grad_deriv)
348 : ! VERY UGLY HACK e_0 has to multiplied by the factor 2
349 6256743 : e_0_br = 2.0_dp*e_0_br
350 : END IF
351 : END IF
352 :
353 : ! ** Now we calculate the pbe cutoff part
354 : ! ** Attention we need to scale rho, ndrho first
355 6400000 : my_rho = my_rho*2.0_dp
356 6400000 : my_ndrho = my_ndrho*2.0_dp
357 :
358 : ! ** Do some precalculation in order to catch the correct branch afterwards
359 6400000 : sscale = 1.0_dp
360 6400000 : t1 = pi**2
361 6400000 : t2 = t1*my_rho
362 6400000 : t3 = t2**(0.1e1_dp/0.3e1_dp)
363 6400000 : t4 = 0.1e1_dp/t3
364 6400000 : t6 = my_ndrho*t4
365 6400000 : t7 = 0.1e1_dp/my_rho
366 6400000 : t8 = t7*sscale
367 6400000 : ss = 0.3466806371753173524216762e0_dp*t6*t8
368 6400000 : IF (ss > scutoff) THEN
369 3573347 : ss2 = ss*ss
370 3573347 : sscale = (smax*ss2 - sconst)/(ss2*ss)
371 : END IF
372 6400000 : e_0_pbe = 0.0_dp
373 6400000 : e_rho_pbe = 0.0_dp
374 6400000 : e_ndrho_pbe = 0.0_dp
375 6400000 : IF (ss*sscale > gcutoff) THEN
376 : CALL xpbe_hole_t_c_lr_lda_calc_1(e_0_pbe, e_rho_pbe, e_ndrho_pbe, &
377 : my_rho, &
378 6399421 : my_ndrho, sscale, sx, R, grad_deriv)
379 : ELSE
380 : CALL xpbe_hole_t_c_lr_lda_calc_2(e_0_pbe, e_rho_pbe, e_ndrho_pbe, &
381 : my_rho, &
382 579 : my_ndrho, sscale, sx, R, grad_deriv)
383 : END IF
384 :
385 : ! ** Finally we get the LDA part
386 :
387 6400000 : e_0_lda = 0.0_dp
388 6400000 : e_rho_lda = 0.0_dp
389 : CALL xlda_hole_t_c_lr_lda_calc_0(grad_deriv, my_rho, e_0_lda, e_rho_lda, &
390 6400000 : sx, R)
391 :
392 6400000 : Fx = e_0_br/e_0_lda
393 :
394 6400000 : Fermi = alpha/(EXP((Fx - mu)/N) + 1.0_dp)
395 :
396 6400000 : dFermi_drho = -Fermi**2/alpha/N*(e_rho_br/e_0_lda - e_0_br*e_rho_lda/e_0_lda**2)*EXP((Fx - mu)/N)
397 6400000 : dFermi_dndrho = -Fermi**2/alpha/N*(e_ndrho_br/e_0_lda)*EXP((Fx - mu)/N)
398 6400000 : dFermi_dtau = -Fermi**2/alpha/N*(e_tau_br/e_0_lda)*EXP((Fx - mu)/N)
399 6400000 : dFermi_dlaplace_rho = -Fermi**2/alpha/N*(e_laplace_rho_br/e_0_lda)*EXP((Fx - mu)/N)
400 :
401 6400000 : e_0(ip) = e_0(ip) + (Fermi*e_0_pbe + (1.0_dp - Fermi)*e_0_br)*sx
402 :
403 6400000 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
404 :
405 : e_rho(ip) = e_rho(ip) + (Fermi*e_rho_pbe + dFermi_drho*e_0_pbe + &
406 3264000 : (1.0_dp - Fermi)*e_rho_br - dFermi_drho*e_0_br)*sx
407 :
408 : e_ndrho(ip) = e_ndrho(ip) + (Fermi*e_ndrho_pbe + dFermi_dndrho*e_0_pbe + &
409 3264000 : (1.0_dp - Fermi)*e_ndrho_br - dFermi_dndrho*e_0_br)*sx
410 :
411 : e_tau(ip) = e_tau(ip) + (dFermi_dtau*e_0_pbe + &
412 3264000 : (1.0_dp - Fermi)*e_tau_br - dFermi_dtau*e_0_br)*sx
413 :
414 : e_laplace_rho(ip) = e_laplace_rho(ip) + (dFermi_dlaplace_rho*e_0_pbe + &
415 3264000 : (1.0_dp - Fermi)*e_laplace_rho_br - dFermi_dlaplace_rho*e_0_br)*sx
416 : END IF
417 :
418 : END IF
419 : END DO
420 :
421 : !$OMP END DO
422 :
423 200 : END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lda_calc
424 :
425 : ! **************************************************************************************************
426 : !> \brief Intermediate routine that gets grids, derivatives and some params
427 : !> \param rho_set the density where you want to evaluate the functional
428 : !> \param deriv_set place where to store the functional derivatives (they are
429 : !> added to the derivatives)
430 : !> \param grad_deriv degree of the derivative that should be evaluated,
431 : !> if positive all the derivatives up to the given degree are evaluated,
432 : !> if negative only the given degree is calculated
433 : !> \param params parameters for functional
434 : !> \author mguidon (01.2009)
435 : ! **************************************************************************************************
436 290 : SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_eval(rho_set, deriv_set, grad_deriv, params)
437 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
438 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
439 : INTEGER, INTENT(in) :: grad_deriv
440 : TYPE(section_vals_type), POINTER :: params
441 :
442 : CHARACTER(len=*), PARAMETER :: routineN = 'xbr_pbe_lda_hole_tc_lr_lsd_eval'
443 :
444 : INTEGER :: handle, npoints
445 : INTEGER, DIMENSION(2, 3) :: bo
446 : REAL(dp) :: gamma, R, sx
447 : REAL(kind=dp) :: epsilon_rho
448 58 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_laplace_rhoa, &
449 58 : e_laplace_rhob, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, e_tau_a, e_tau_b, laplace_rhoa, &
450 58 : laplace_rhob, norm_drhoa, norm_drhob, rhoa, rhob, tau_a, tau_b
451 : TYPE(xc_derivative_type), POINTER :: deriv
452 :
453 58 : CALL timeset(routineN, handle)
454 :
455 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
456 : norm_drhob=norm_drhob, tau_a=tau_a, tau_b=tau_b, laplace_rhoa=laplace_rhoa, &
457 : laplace_rhob=laplace_rhob, local_bounds=bo, &
458 58 : rho_cutoff=epsilon_rho)
459 58 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
460 :
461 58 : dummy => rhoa
462 :
463 58 : e_0 => dummy
464 58 : e_rhoa => dummy
465 58 : e_rhob => dummy
466 58 : e_ndrhoa => dummy
467 58 : e_ndrhob => dummy
468 58 : e_tau_a => dummy
469 58 : e_tau_b => dummy
470 58 : e_laplace_rhoa => dummy
471 58 : e_laplace_rhob => dummy
472 :
473 58 : IF (grad_deriv >= 0) THEN
474 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
475 58 : allocate_deriv=.TRUE.)
476 58 : CALL xc_derivative_get(deriv, deriv_data=e_0)
477 : END IF
478 58 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
479 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
480 30 : allocate_deriv=.TRUE.)
481 30 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
482 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
483 30 : allocate_deriv=.TRUE.)
484 30 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
485 :
486 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
487 30 : allocate_deriv=.TRUE.)
488 30 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
489 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
490 30 : allocate_deriv=.TRUE.)
491 30 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
492 :
493 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], &
494 30 : allocate_deriv=.TRUE.)
495 30 : CALL xc_derivative_get(deriv, deriv_data=e_tau_a)
496 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], &
497 30 : allocate_deriv=.TRUE.)
498 30 : CALL xc_derivative_get(deriv, deriv_data=e_tau_b)
499 :
500 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhoa], &
501 30 : allocate_deriv=.TRUE.)
502 30 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhoa)
503 : deriv => xc_dset_get_derivative(deriv_set, [deriv_laplace_rhob], &
504 30 : allocate_deriv=.TRUE.)
505 30 : CALL xc_derivative_get(deriv, deriv_data=e_laplace_rhob)
506 : END IF
507 58 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
508 0 : CPABORT("derivatives bigger than 1 not implemented")
509 : END IF
510 :
511 58 : CALL section_vals_val_get(params, "scale_x", r_val=sx)
512 58 : CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
513 58 : CALL section_vals_val_get(params, "GAMMA", r_val=gamma)
514 :
515 58 : IF (R == 0.0_dp) THEN
516 0 : CPABORT("Cutoff_Radius 0.0 not implemented")
517 : END IF
518 :
519 : !$OMP PARALLEL DEFAULT(NONE) &
520 : !$OMP SHARED(rhoa, norm_drhoa, laplace_rhoa, tau_a, e_0) &
521 : !$OMP SHARED(e_rhoa, e_ndrhoa, e_tau_a, e_laplace_rhoa) &
522 : !$OMP SHARED(grad_deriv, npoints, epsilon_rho) &
523 : !$OMP SHARED(sx, r, gamma) &
524 : !$OMP SHARED(rhob, norm_drhob, laplace_rhob, tau_b, e_rhob) &
525 58 : !$OMP SHARED(e_ndrhob, e_tau_b, e_laplace_rhob)
526 :
527 : CALL xbr_pbe_lda_hole_tc_lr_lsd_calc(rho=rhoa, norm_drho=norm_drhoa, &
528 : laplace_rho=laplace_rhoa, tau=tau_a, e_0=e_0, e_rho=e_rhoa, e_ndrho=e_ndrhoa, &
529 : e_tau=e_tau_a, e_laplace_rho=e_laplace_rhoa, grad_deriv=grad_deriv, &
530 : npoints=npoints, epsilon_rho=epsilon_rho, &
531 : sx=sx, R=R, gamma=gamma)
532 :
533 : CALL xbr_pbe_lda_hole_tc_lr_lsd_calc(rho=rhob, norm_drho=norm_drhob, &
534 : laplace_rho=laplace_rhob, tau=tau_b, e_0=e_0, e_rho=e_rhob, e_ndrho=e_ndrhob, &
535 : e_tau=e_tau_b, e_laplace_rho=e_laplace_rhob, grad_deriv=grad_deriv, &
536 : npoints=npoints, epsilon_rho=epsilon_rho, &
537 : sx=sx, R=R, gamma=gamma)
538 :
539 : !$OMP END PARALLEL
540 :
541 58 : CALL timestop(handle)
542 58 : END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_eval
543 : ! **************************************************************************************************
544 : !> \brief Low level routine that calls the three involved holes and puts them
545 : !> together
546 : !> \param rho values on the grid
547 : !> \param norm_drho values on the grid
548 : !> \param laplace_rho values on the grid
549 : !> \param tau values on the grid
550 : !> \param e_0 derivatives on the grid
551 : !> \param e_rho derivatives on the grid
552 : !> \param e_ndrho derivatives on the grid
553 : !> \param e_tau derivatives on the grid
554 : !> \param e_laplace_rho derivatives on the grid
555 : !> \param grad_deriv degree of the derivative that should be evaluated,
556 : !> if positive all the derivatives up to the given degree are evaluated,
557 : !> if negative only the given degree is calculated
558 : !> \param npoints number of gridpoints
559 : !> \param epsilon_rho cutoffs
560 : !> \param sx parameters for functional
561 : !> \param R parameters for functional
562 : !> \param gamma parameters for functional
563 : !> \author mguidon (01.2009)
564 : ! **************************************************************************************************
565 116 : SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_calc(rho, norm_drho, laplace_rho, tau, e_0, e_rho, &
566 116 : e_ndrho, e_tau, e_laplace_rho, grad_deriv, npoints, &
567 : epsilon_rho, sx, R, gamma)
568 :
569 : INTEGER, INTENT(in) :: npoints, grad_deriv
570 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_laplace_rho, e_tau, e_ndrho, e_rho, e_0
571 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: tau, laplace_rho, norm_drho, rho
572 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R, gamma
573 :
574 : INTEGER :: ip
575 : REAL(dp) :: dFermi_dlaplace_rho, dFermi_dndrho, dFermi_drho, dFermi_dtau, e_0_br, e_0_lda, &
576 : e_0_pbe, e_laplace_rho_br, e_ndrho_br, e_ndrho_pbe, e_rho_br, e_rho_lda, e_rho_pbe, &
577 : e_tau_br, Fermi, Fx, my_laplace_rho, my_ndrho, my_rho, my_tau, ss, ss2, sscale, t1, t15, &
578 : t16, t2, t3, t4, t5, t6, t7, t8, t9, yval
579 :
580 : !$OMP DO
581 :
582 : DO ip = 1, npoints
583 5285250 : my_rho = MAX(rho(ip), 0.0_dp)
584 5285250 : IF (my_rho > epsilon_rho) THEN
585 5285248 : my_ndrho = MAX(norm_drho(ip), EPSILON(0.0_dp)*1.e4_dp)
586 5285248 : my_tau = 1.0_dp*MAX(EPSILON(0.0_dp)*1.e4_dp, tau(ip))
587 5285248 : my_laplace_rho = 1.0_dp*laplace_rho(ip)
588 :
589 5285248 : t1 = pi**(0.1e1_dp/0.3e1_dp)
590 5285248 : t2 = t1**2
591 5285248 : t3 = my_rho**(0.1e1_dp/0.3e1_dp)
592 5285248 : t4 = t3**2
593 5285248 : t5 = t4*my_rho
594 5285248 : t8 = my_ndrho**2
595 5285248 : t9 = 0.1e1_dp/my_rho
596 5285248 : t15 = my_laplace_rho/0.6e1_dp - gamma*(2.0_dp*my_tau - t8*t9/0.4e1_dp)/0.3e1_dp
597 5285248 : t16 = 0.1e1_dp/t15
598 5285248 : yval = 0.2e1_dp/0.3e1_dp*t2*t5*t16
599 :
600 5285248 : e_0_br = 0.0_dp
601 5285248 : e_rho_br = 0.0_dp
602 5285248 : e_ndrho_br = 0.0_dp
603 5285248 : e_tau_br = 0.0_dp
604 5285248 : e_laplace_rho_br = 0.0_dp
605 :
606 5285248 : IF (R == 0.0_dp) THEN
607 0 : IF (yval <= 0.0_dp) THEN
608 : CALL x_br_lsd_y_lte_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
609 : e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
610 0 : sx, gamma, grad_deriv)
611 : ELSE
612 : CALL x_br_lsd_y_gt_0(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
613 : e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
614 0 : sx, gamma, grad_deriv)
615 : END IF
616 : ELSE
617 5285248 : IF (yval <= 0.0_dp) THEN
618 : CALL x_br_lsd_y_lte_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
619 : e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
620 115013 : sx, R, gamma, grad_deriv)
621 : ELSE
622 : CALL x_br_lsd_y_gt_0_cutoff(my_rho, my_ndrho, my_tau, my_laplace_rho, e_0_br, &
623 : e_rho_br, e_ndrho_br, e_tau_br, e_laplace_rho_br, &
624 5170235 : sx, R, gamma, grad_deriv)
625 : END IF
626 : END IF
627 :
628 : ! ** Now we calculate the pbe cutoff part
629 : ! ** Attention we need to scale rho, ndrho first
630 5285248 : my_rho = my_rho*2.0_dp
631 5285248 : my_ndrho = my_ndrho*2.0_dp
632 :
633 : ! ** Do some precalculation in order to catch the correct branch afterwards
634 5285248 : sscale = 1.0_dp
635 5285248 : t1 = pi**2
636 5285248 : t2 = t1*my_rho
637 5285248 : t3 = t2**(0.1e1_dp/0.3e1_dp)
638 5285248 : t4 = 0.1e1_dp/t3
639 5285248 : t6 = my_ndrho*t4
640 5285248 : t7 = 0.1e1_dp/my_rho
641 5285248 : t8 = t7*sscale
642 5285248 : ss = 0.3466806371753173524216762e0_dp*t6*t8
643 5285248 : IF (ss > scutoff) THEN
644 1684625 : ss2 = ss*ss
645 1684625 : sscale = (smax*ss2 - sconst)/(ss2*ss)
646 : END IF
647 5285248 : e_0_pbe = 0.0_dp
648 5285248 : e_rho_pbe = 0.0_dp
649 5285248 : e_ndrho_pbe = 0.0_dp
650 5285248 : IF (ss*sscale > gcutoff) THEN
651 : CALL xpbe_hole_t_c_lr_lda_calc_1(e_0_pbe, e_rho_pbe, e_ndrho_pbe, &
652 : my_rho, &
653 5285085 : my_ndrho, sscale, sx, R, grad_deriv)
654 : ELSE
655 : CALL xpbe_hole_t_c_lr_lda_calc_2(e_0_pbe, e_rho_pbe, e_ndrho_pbe, &
656 : my_rho, &
657 163 : my_ndrho, sscale, sx, R, grad_deriv)
658 : END IF
659 :
660 5285248 : e_0_pbe = 0.5_dp*e_0_pbe
661 :
662 : ! ** Finally we get the LDA part
663 :
664 5285248 : e_0_lda = 0.0_dp
665 5285248 : e_rho_lda = 0.0_dp
666 : CALL xlda_hole_t_c_lr_lda_calc_0(grad_deriv, my_rho, e_0_lda, e_rho_lda, &
667 5285248 : sx, R)
668 5285248 : e_0_lda = 0.5_dp*e_0_lda
669 :
670 5285248 : Fx = e_0_br/e_0_lda
671 :
672 5285248 : Fermi = alpha/(EXP((Fx - mu)/N) + 1.0_dp)
673 :
674 5285248 : dFermi_drho = -Fermi**2/alpha/N*(e_rho_br/e_0_lda - e_0_br*e_rho_lda/e_0_lda**2)*EXP((Fx - mu)/N)
675 5285248 : dFermi_dndrho = -Fermi**2/alpha/N*(e_ndrho_br/e_0_lda)*EXP((Fx - mu)/N)
676 5285248 : dFermi_dtau = -Fermi**2/alpha/N*(e_tau_br/e_0_lda)*EXP((Fx - mu)/N)
677 5285248 : dFermi_dlaplace_rho = -Fermi**2/alpha/N*(e_laplace_rho_br/e_0_lda)*EXP((Fx - mu)/N)
678 :
679 5285248 : e_0(ip) = e_0(ip) + (Fermi*e_0_pbe + (1.0_dp - Fermi)*e_0_br)*sx
680 :
681 5285248 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
682 :
683 : e_rho(ip) = e_rho(ip) + (Fermi*e_rho_pbe + dFermi_drho*e_0_pbe + &
684 2733750 : (1.0_dp - Fermi)*e_rho_br - dFermi_drho*e_0_br)*sx
685 :
686 : e_ndrho(ip) = e_ndrho(ip) + (Fermi*e_ndrho_pbe + dFermi_dndrho*e_0_pbe + &
687 2733750 : (1.0_dp - Fermi)*e_ndrho_br - dFermi_dndrho*e_0_br)*sx
688 :
689 : e_tau(ip) = e_tau(ip) + (dFermi_dtau*e_0_pbe + &
690 2733750 : (1.0_dp - Fermi)*e_tau_br - dFermi_dtau*e_0_br)*sx
691 :
692 : e_laplace_rho(ip) = e_laplace_rho(ip) + (dFermi_dlaplace_rho*e_0_pbe + &
693 2733750 : (1.0_dp - Fermi)*e_laplace_rho_br - dFermi_dlaplace_rho*e_0_br)*sx
694 : END IF
695 :
696 : END IF
697 : END DO
698 :
699 : !$OMP END DO
700 :
701 116 : END SUBROUTINE xbr_pbe_lda_hole_tc_lr_lsd_calc
702 :
703 : END MODULE xc_xbr_pbe_lda_hole_t_c_lr
|
