Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief calculates the Becke 88 exchange functional
10 : !> \par History
11 : !> 11.2003 created [fawzi]
12 : !> \author fawzi
13 : ! **************************************************************************************************
14 : MODULE xc_xbecke88
15 : USE bibliography, ONLY: Becke1988,&
16 : cite_reference
17 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
18 : USE input_section_types, ONLY: section_vals_type,&
19 : section_vals_val_get
20 : USE kinds, ONLY: dp
21 : USE mathconstants, ONLY: pi
22 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
23 : deriv_norm_drhoa,&
24 : deriv_norm_drhob,&
25 : deriv_rho,&
26 : deriv_rhoa,&
27 : deriv_rhob
28 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
29 : xc_dset_get_derivative
30 : USE xc_derivative_types, ONLY: xc_derivative_get,&
31 : xc_derivative_type
32 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
33 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
34 : xc_rho_set_type
35 : #include "../base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 : PRIVATE
39 :
40 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
41 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xbecke88'
42 : REAL(kind=dp), PARAMETER :: beta = 0.0042_dp
43 :
44 : PUBLIC :: xb88_lda_info, xb88_lsd_info, xb88_lda_eval, xb88_lsd_eval
45 : CONTAINS
46 :
47 : ! **************************************************************************************************
48 : !> \brief return various information on the functional
49 : !> \param reference string with the reference of the actual functional
50 : !> \param shortform string with the shortform of the functional name
51 : !> \param needs the components needed by this functional are set to
52 : !> true (does not set the unneeded components to false)
53 : !> \param max_deriv ...
54 : !> \par History
55 : !> 11.2003 created [fawzi]
56 : !> \author fawzi
57 : ! **************************************************************************************************
58 9702 : SUBROUTINE xb88_lda_info(reference, shortform, needs, max_deriv)
59 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
60 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
61 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
62 :
63 9702 : IF (PRESENT(reference)) THEN
64 91 : reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LDA version}"
65 : END IF
66 9702 : IF (PRESENT(shortform)) THEN
67 91 : shortform = "Becke 1988 Exchange Functional (LDA)"
68 : END IF
69 9702 : IF (PRESENT(needs)) THEN
70 9611 : needs%rho = .TRUE.
71 9611 : needs%rho_1_3 = .TRUE.
72 9611 : needs%norm_drho = .TRUE.
73 : END IF
74 9702 : IF (PRESENT(max_deriv)) max_deriv = 3
75 :
76 9702 : END SUBROUTINE xb88_lda_info
77 :
78 : ! **************************************************************************************************
79 : !> \brief return various information on the functional
80 : !> \param reference string with the reference of the actual functional
81 : !> \param shortform string with the shortform of the functional name
82 : !> \param needs the components needed by this functional are set to
83 : !> true (does not set the unneeded components to false)
84 : !> \param max_deriv ...
85 : !> \par History
86 : !> 11.2003 created [fawzi]
87 : !> \author fawzi
88 : ! **************************************************************************************************
89 2839 : SUBROUTINE xb88_lsd_info(reference, shortform, needs, max_deriv)
90 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
91 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
92 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
93 :
94 2839 : IF (PRESENT(reference)) THEN
95 41 : reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LSD version}"
96 : END IF
97 2839 : IF (PRESENT(shortform)) THEN
98 41 : shortform = "Becke 1988 Exchange Functional (LSD)"
99 : END IF
100 2839 : IF (PRESENT(needs)) THEN
101 2798 : needs%rho_spin = .TRUE.
102 2798 : needs%rho_spin_1_3 = .TRUE.
103 2798 : needs%norm_drho_spin = .TRUE.
104 : END IF
105 2839 : IF (PRESENT(max_deriv)) max_deriv = 3
106 :
107 2839 : END SUBROUTINE xb88_lsd_info
108 :
109 : ! **************************************************************************************************
110 : !> \brief evaluates the becke 88 exchange functional for lda
111 : !> \param rho_set the density where you want to evaluate the functional
112 : !> \param deriv_set place where to store the functional derivatives (they are
113 : !> added to the derivatives)
114 : !> \param grad_deriv degree of the derivative that should be evaluated,
115 : !> if positive all the derivatives up to the given degree are evaluated,
116 : !> if negative only the given degree is calculated
117 : !> \param xb88_params input parameters (scaling)
118 : !> \par History
119 : !> 11.2003 created [fawzi]
120 : !> 01.2007 added scaling [Manuel Guidon]
121 : !> \author fawzi
122 : ! **************************************************************************************************
123 27171 : SUBROUTINE xb88_lda_eval(rho_set, deriv_set, grad_deriv, xb88_params)
124 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
125 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
126 : INTEGER, INTENT(in) :: grad_deriv
127 : TYPE(section_vals_type), POINTER :: xb88_params
128 :
129 : CHARACTER(len=*), PARAMETER :: routineN = 'xb88_lda_eval'
130 :
131 : INTEGER :: handle, npoints
132 : INTEGER, DIMENSION(2, 3) :: bo
133 : REAL(kind=dp) :: epsilon_rho, sx
134 9057 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
135 9057 : e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
136 9057 : e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho, rho_1_3
137 : TYPE(xc_derivative_type), POINTER :: deriv
138 :
139 9057 : CALL timeset(routineN, handle)
140 :
141 9057 : CALL section_vals_val_get(xb88_params, "scale_x", r_val=sx)
142 :
143 9057 : CALL cite_reference(Becke1988)
144 :
145 : CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
146 9057 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
147 9057 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
148 :
149 9057 : dummy => rho
150 :
151 9057 : e_0 => dummy
152 9057 : e_rho => dummy
153 9057 : e_ndrho => dummy
154 9057 : e_rho_rho => dummy
155 9057 : e_ndrho_rho => dummy
156 9057 : e_ndrho_ndrho => dummy
157 9057 : e_rho_rho_rho => dummy
158 9057 : e_ndrho_rho_rho => dummy
159 9057 : e_ndrho_ndrho_rho => dummy
160 9057 : e_ndrho_ndrho_ndrho => dummy
161 :
162 9057 : IF (grad_deriv >= 0) THEN
163 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
164 9057 : allocate_deriv=.TRUE.)
165 9057 : CALL xc_derivative_get(deriv, deriv_data=e_0)
166 : END IF
167 9057 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
168 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
169 8815 : allocate_deriv=.TRUE.)
170 8815 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
171 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
172 8815 : allocate_deriv=.TRUE.)
173 8815 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
174 : END IF
175 9057 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
176 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
177 628 : allocate_deriv=.TRUE.)
178 628 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
179 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
180 628 : allocate_deriv=.TRUE.)
181 628 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
182 : deriv => xc_dset_get_derivative(deriv_set, &
183 628 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
184 628 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
185 : END IF
186 9057 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
187 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
188 0 : allocate_deriv=.TRUE.)
189 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
190 : deriv => xc_dset_get_derivative(deriv_set, &
191 0 : [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.TRUE.)
192 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
193 : deriv => xc_dset_get_derivative(deriv_set, &
194 0 : [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.TRUE.)
195 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
196 : deriv => xc_dset_get_derivative(deriv_set, &
197 0 : [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
198 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
199 : END IF
200 9057 : IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
201 0 : CPABORT("derivatives bigger than 3 not implemented")
202 : END IF
203 :
204 : !$OMP PARALLEL DEFAULT(NONE) &
205 : !$OMP SHARED(rho, rho_1_3, norm_drho, e_0, e_rho) &
206 : !$OMP SHARED(e_ndrho, e_rho_rho, e_ndrho_rho) &
207 : !$OMP SHARED(e_ndrho_ndrho, e_rho_rho_rho) &
208 : !$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
209 : !$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
210 9057 : !$OMP SHARED(epsilon_rho, sx)
211 :
212 : CALL xb88_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
213 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
214 : e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
215 : e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
216 : e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
217 : e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, grad_deriv=grad_deriv, &
218 : npoints=npoints, epsilon_rho=epsilon_rho, sx=sx)
219 :
220 : !$OMP END PARALLEL
221 :
222 9057 : CALL timestop(handle)
223 9057 : END SUBROUTINE xb88_lda_eval
224 :
225 : ! **************************************************************************************************
226 : !> \brief evaluates the becke 88 exchange functional for lda
227 : !> \param rho the density where you want to evaluate the functional
228 : !> \param rho_1_3 ...
229 : !> \param norm_drho ...
230 : !> \param e_0 ...
231 : !> \param e_rho ...
232 : !> \param e_ndrho ...
233 : !> \param e_rho_rho ...
234 : !> \param e_ndrho_rho ...
235 : !> \param e_ndrho_ndrho ...
236 : !> \param e_rho_rho_rho ...
237 : !> \param e_ndrho_rho_rho ...
238 : !> \param e_ndrho_ndrho_rho ...
239 : !> \param e_ndrho_ndrho_ndrho ...
240 : !> \param grad_deriv degree of the derivative that should be evaluated,
241 : !> if positive all the derivatives up to the given degree are evaluated,
242 : !> if negative only the given degree is calculated
243 : !> \param npoints ...
244 : !> \param epsilon_rho ...
245 : !> \param sx scaling-parameter for exchange
246 : !> \par History
247 : !> 11.2003 created [fawzi]
248 : !> 01.2007 added scaling [Manuel Guidon]
249 : !> \author fawzi
250 : ! **************************************************************************************************
251 9057 : SUBROUTINE xb88_lda_calc(rho, rho_1_3, norm_drho, &
252 9057 : e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
253 9057 : e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
254 9057 : e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, sx)
255 : INTEGER, INTENT(in) :: npoints, grad_deriv
256 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_ndrho, &
257 : e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, &
258 : e_ndrho, e_rho, e_0
259 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho_1_3, rho
260 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
261 :
262 : INTEGER :: ii
263 : REAL(kind=dp) :: c, epsilon_rho43, my_rho, my_rho_1_3, t1, t10, t100, t104, t11, t12, t126, &
264 : t13, t159, t164, t17, t170, t176, t18, t189, t2, t21, t23, t29, t3, t31, t33, t37, t39, &
265 : t40, t43, t44, t49, t5, t51, t54, t6, t64, t67, t7, t72, t74, t75, t79, t83, t86, t90, &
266 : t91, t98, t99, x
267 :
268 9057 : c = 1.5_dp*(0.75_dp/Pi)**(1.0_dp/3.0_dp)
269 9057 : t1 = 2**(0.1e1_dp/0.3e1_dp)
270 9057 : t2 = t1**2
271 9057 : t5 = beta*t2
272 9057 : t7 = beta*t1
273 9057 : epsilon_rho43 = epsilon_rho**(4./3.)
274 :
275 242 : SELECT CASE (grad_deriv)
276 : CASE (0)
277 :
278 8187 : !$OMP DO
279 :
280 : DO ii = 1, npoints
281 :
282 28923624 : my_rho = rho(ii)
283 28923624 : IF (my_rho > epsilon_rho) THEN
284 23036823 : my_rho_1_3 = rho_1_3(ii)
285 23036823 : t3 = my_rho_1_3*my_rho
286 23036823 : x = norm_drho(ii)/MAX(t3, epsilon_rho43)
287 23036823 : t6 = x**2
288 23036823 : t10 = t2*t6 + 0.1e1_dp
289 23036823 : t11 = SQRT(t10)
290 23036823 : t12 = t1*x + t11
291 23036823 : t13 = LOG(t12)
292 23036823 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
293 23036823 : t18 = 0.1e1_dp/t17
294 23036823 : t21 = -c - t5*t6*t18
295 :
296 23036823 : e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
297 : END IF
298 : END DO
299 :
300 : !$OMP END DO
301 :
302 : CASE (1)
303 :
304 0 : !$OMP DO
305 :
306 : DO ii = 1, npoints
307 :
308 323835149 : my_rho = rho(ii)
309 323835149 : IF (my_rho > epsilon_rho) THEN
310 310495634 : my_rho_1_3 = rho_1_3(ii)
311 310495634 : t3 = my_rho_1_3*my_rho
312 310495634 : x = norm_drho(ii)/t3
313 310495634 : t6 = x**2
314 310495634 : t10 = t2*t6 + 0.1e1_dp
315 310495634 : t11 = SQRT(t10)
316 310495634 : t12 = t1*x + t11
317 310495634 : t13 = LOG(t12)
318 310495634 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
319 310495634 : t18 = 0.1e1_dp/t17
320 310495634 : t21 = -c - t5*t6*t18
321 :
322 310495634 : e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
323 :
324 310495634 : t23 = t2*my_rho_1_3
325 310495634 : t29 = 0.1e1_dp/t11
326 310495634 : t31 = t1 + t2*x*t29
327 310495634 : t33 = 0.1e1_dp/t12
328 310495634 : t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
329 310495634 : t39 = t17**2
330 310495634 : t40 = 0.1e1_dp/t39
331 310495634 : t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
332 310495634 : t44 = t43*x
333 :
334 : e_rho(ii) = e_rho(ii) &
335 : - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
336 310495634 : t23*t21)*sx
337 : e_ndrho(ii) = e_ndrho(ii) &
338 310495634 : + (t2*t43/0.2e1_dp)*sx
339 : END IF
340 :
341 : END DO
342 :
343 : !$OMP END DO
344 :
345 : CASE (-1)
346 :
347 628 : !$OMP DO
348 :
349 : DO ii = 1, npoints
350 :
351 0 : my_rho = rho(ii)
352 0 : IF (my_rho > epsilon_rho) THEN
353 0 : my_rho_1_3 = rho_1_3(ii)
354 0 : t3 = my_rho_1_3*my_rho
355 0 : x = norm_drho(ii)/t3
356 0 : t6 = x**2
357 0 : t10 = t2*t6 + 0.1e1_dp
358 0 : t11 = SQRT(t10)
359 0 : t12 = t1*x + t11
360 0 : t13 = LOG(t12)
361 0 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
362 0 : t18 = 0.1e1_dp/t17
363 0 : t21 = -c - t5*t6*t18
364 :
365 0 : t23 = t2*my_rho_1_3
366 0 : t29 = 0.1e1_dp/t11
367 0 : t31 = t1 + t2*x*t29
368 0 : t33 = 0.1e1_dp/t12
369 0 : t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
370 0 : t39 = t17**2
371 0 : t40 = 0.1e1_dp/t39
372 0 : t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
373 0 : t44 = t43*x
374 :
375 : e_rho(ii) = e_rho(ii) &
376 : - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
377 0 : t23*t21)*sx
378 : e_ndrho(ii) = e_ndrho(ii) &
379 0 : + (t2*t43/0.2e1_dp)*sx
380 : END IF
381 : END DO
382 :
383 : !$OMP END DO
384 :
385 : CASE (2)
386 :
387 0 : !$OMP DO
388 :
389 : DO ii = 1, npoints
390 :
391 13835897 : my_rho = rho(ii)
392 13835897 : IF (my_rho > epsilon_rho) THEN
393 13647095 : my_rho_1_3 = rho_1_3(ii)
394 13647095 : t3 = my_rho_1_3*my_rho
395 13647095 : x = norm_drho(ii)/t3
396 13647095 : t6 = x**2
397 13647095 : t10 = t2*t6 + 0.1e1_dp
398 13647095 : t11 = SQRT(t10)
399 13647095 : t12 = t1*x + t11
400 13647095 : t13 = LOG(t12)
401 13647095 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
402 13647095 : t18 = 0.1e1_dp/t17
403 13647095 : t21 = -c - t5*t6*t18
404 :
405 13647095 : e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
406 :
407 13647095 : t23 = t2*my_rho_1_3
408 13647095 : t29 = 0.1e1_dp/t11
409 13647095 : t31 = t1 + t2*x*t29
410 13647095 : t33 = 0.1e1_dp/t12
411 13647095 : t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
412 13647095 : t39 = t17**2
413 13647095 : t40 = 0.1e1_dp/t39
414 13647095 : t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
415 13647095 : t44 = t43*x
416 :
417 : e_rho(ii) = e_rho(ii) &
418 : - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
419 13647095 : t23*t21)*sx
420 : e_ndrho(ii) = e_ndrho(ii) &
421 13647095 : + (t2*t43/0.2e1_dp)*sx
422 :
423 13647095 : t49 = my_rho_1_3**2
424 13647095 : t51 = t2/t49
425 13647095 : t54 = x*t40
426 13647095 : t64 = 0.1e1_dp/t11/t10
427 13647095 : t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
428 13647095 : t72 = t31**2
429 13647095 : t74 = t12**2
430 13647095 : t75 = 0.1e1_dp/t74
431 : t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
432 13647095 : - 0.6e1_dp*t7*x*t72*t75
433 13647095 : t83 = t37**2
434 13647095 : t86 = 0.1e1_dp/t39/t17
435 : t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
436 13647095 : t79*t40 - 0.2e1_dp*t5*t6*t83*t86
437 13647095 : t91 = t90*t6
438 13647095 : t98 = 0.1e1_dp/my_rho
439 13647095 : t99 = t2*t98
440 13647095 : t100 = t90*x
441 13647095 : t104 = 0.1e1_dp/t3
442 :
443 : e_rho_rho(ii) = e_rho_rho(ii) &
444 : + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
445 13647095 : t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
446 : e_ndrho_rho(ii) = e_ndrho_rho(ii) &
447 13647095 : - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
448 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
449 13647095 : + (t2*t90*t104/0.2e1_dp)*sx
450 : END IF
451 : END DO
452 :
453 : !$OMP END DO
454 :
455 : CASE (-2)
456 :
457 0 : !$OMP DO
458 :
459 : DO ii = 1, npoints
460 :
461 0 : my_rho = rho(ii)
462 0 : IF (my_rho > epsilon_rho) THEN
463 0 : my_rho_1_3 = rho_1_3(ii)
464 0 : t3 = my_rho_1_3*my_rho
465 0 : x = norm_drho(ii)/t3
466 0 : t6 = x**2
467 0 : t10 = t2*t6 + 0.1e1_dp
468 0 : t11 = SQRT(t10)
469 0 : t12 = t1*x + t11
470 0 : t13 = LOG(t12)
471 0 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
472 0 : t18 = 0.1e1_dp/t17
473 0 : t21 = -c - t5*t6*t18
474 :
475 0 : t23 = t2*my_rho_1_3
476 0 : t29 = 0.1e1_dp/t11
477 0 : t31 = t1 + t2*x*t29
478 0 : t33 = 0.1e1_dp/t12
479 0 : t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
480 0 : t39 = t17**2
481 0 : t40 = 0.1e1_dp/t39
482 0 : t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
483 0 : t44 = t43*x
484 :
485 0 : t49 = my_rho_1_3**2
486 0 : t51 = t2/t49
487 0 : t54 = x*t40
488 0 : t64 = 0.1e1_dp/t11/t10
489 0 : t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
490 0 : t72 = t31**2
491 0 : t74 = t12**2
492 0 : t75 = 0.1e1_dp/t74
493 : t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
494 0 : - 0.6e1_dp*t7*x*t72*t75
495 0 : t83 = t37**2
496 0 : t86 = 0.1e1_dp/t39/t17
497 : t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
498 0 : t79*t40 - 0.2e1_dp*t5*t6*t83*t86
499 0 : t91 = t90*t6
500 0 : t98 = 0.1e1_dp/my_rho
501 0 : t99 = t2*t98
502 0 : t100 = t90*x
503 0 : t104 = 0.1e1_dp/t3
504 :
505 : e_rho_rho(ii) = e_rho_rho(ii) &
506 : + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
507 0 : t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
508 : e_ndrho_rho(ii) = e_ndrho_rho(ii) &
509 0 : - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
510 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
511 0 : + (t2*t90*t104/0.2e1_dp)*sx
512 : END IF
513 : END DO
514 :
515 : !$OMP END DO
516 :
517 : CASE default
518 :
519 9057 : !$OMP DO
520 :
521 : DO ii = 1, npoints
522 :
523 0 : my_rho = rho(ii)
524 0 : IF (my_rho > epsilon_rho) THEN
525 0 : my_rho_1_3 = rho_1_3(ii)
526 0 : t3 = my_rho_1_3*my_rho
527 0 : x = norm_drho(ii)/t3
528 0 : t6 = x**2
529 0 : t10 = t2*t6 + 0.1e1_dp
530 0 : t11 = SQRT(t10)
531 0 : t12 = t1*x + t11
532 0 : t13 = LOG(t12)
533 0 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
534 0 : t18 = 0.1e1_dp/t17
535 0 : t21 = -c - t5*t6*t18
536 :
537 0 : IF (grad_deriv >= 0) THEN
538 0 : e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
539 : END IF
540 :
541 0 : t23 = t2*my_rho_1_3
542 0 : t29 = 0.1e1_dp/t11
543 0 : t31 = t1 + t2*x*t29
544 0 : t33 = 0.1e1_dp/t12
545 0 : t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
546 0 : t39 = t17**2
547 0 : t40 = 0.1e1_dp/t39
548 0 : t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
549 0 : t44 = t43*x
550 :
551 0 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
552 : e_rho(ii) = e_rho(ii) &
553 : - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
554 0 : t23*t21)*sx
555 : e_ndrho(ii) = e_ndrho(ii) &
556 0 : + (t2*t43/0.2e1_dp)*sx
557 : END IF
558 :
559 0 : t49 = my_rho_1_3**2
560 0 : t51 = t2/t49
561 0 : t54 = x*t40
562 0 : t64 = 0.1e1_dp/t11/t10
563 0 : t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
564 0 : t72 = t31**2
565 0 : t74 = t12**2
566 0 : t75 = 0.1e1_dp/t74
567 : t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
568 0 : - 0.6e1_dp*t7*x*t72*t75
569 0 : t83 = t37**2
570 0 : t86 = 0.1e1_dp/t39/t17
571 : t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
572 0 : t79*t40 - 0.2e1_dp*t5*t6*t83*t86
573 0 : t91 = t90*t6
574 0 : t98 = 0.1e1_dp/my_rho
575 0 : t99 = t2*t98
576 0 : t100 = t90*x
577 0 : t104 = 0.1e1_dp/t3
578 :
579 0 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
580 : e_rho_rho(ii) = e_rho_rho(ii) &
581 : + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
582 0 : t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
583 : e_ndrho_rho(ii) = e_ndrho_rho(ii) &
584 0 : - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
585 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
586 0 : + (t2*t90*t104/0.2e1_dp)*sx
587 : END IF
588 :
589 0 : t126 = t10**2
590 0 : t159 = t39**2
591 : t164 = 0.6e1_dp*t5*t40*t37 - 0.12e2_dp*t5*x*t86*t83 &
592 : + 0.6e1_dp*t5*t54*t79 + t5*t6*(0.18e2_dp*t7*t67 &
593 : *t33 - 0.18e2_dp*t7*t72*t75 + 0.6e1_dp*t7*x* &
594 : (-0.6e1_dp*t1*t64*x + 0.12e2_dp*t6*x/t11/t126) &
595 : *t33 - 0.18e2_dp*t7*x*t67*t75*t31 + 0.12e2_dp*t7 &
596 : *x*t72*t31/t74/t12)*t40 - 0.6e1_dp*t5*t6*t79 &
597 0 : *t86*t37 + 0.6e1_dp*t5*t6*t83*t37/t159
598 : t170 = 0.8e1_dp/0.9e1_dp*t51*t164*t6 + 0.14e2_dp/0.9e1_dp &
599 0 : *t51*t100
600 0 : t176 = t2/t49/my_rho
601 0 : t189 = my_rho**2
602 :
603 0 : IF (grad_deriv == -3 .OR. grad_deriv >= 3) THEN
604 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) &
605 : - (0.4e1_dp/0.3e1_dp*t170*x*t98 + 0.16e2_dp/0.27e2_dp &
606 : *t176*t91 - 0.4e1_dp/0.27e2_dp*t176*t44 + &
607 0 : 0.4e1_dp/0.27e2_dp*t176*t21)*sx
608 : e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) &
609 0 : + (t170*t104)*sx
610 : e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) &
611 : + ((-0.2e1_dp/0.3e1_dp*t99*t164*x - &
612 0 : 0.2e1_dp/0.3e1_dp*t99*t90)*t104)*sx
613 : e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) &
614 0 : + (t2*t164/t49/t189/0.2e1_dp)*sx
615 : END IF
616 : END IF
617 : END DO
618 :
619 : !$OMP END DO
620 :
621 : END SELECT
622 :
623 9057 : END SUBROUTINE xb88_lda_calc
624 :
625 : ! **************************************************************************************************
626 : !> \brief evaluates the becke 88 exchange functional for lsd
627 : !> \param rho_set ...
628 : !> \param deriv_set ...
629 : !> \param grad_deriv ...
630 : !> \param xb88_params ...
631 : !> \par History
632 : !> 11.2003 created [fawzi]
633 : !> 01.2007 added scaling [Manuel Guidon]
634 : !> \author fawzi
635 : ! **************************************************************************************************
636 12368 : SUBROUTINE xb88_lsd_eval(rho_set, deriv_set, grad_deriv, xb88_params)
637 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
638 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
639 : INTEGER, INTENT(in) :: grad_deriv
640 : TYPE(section_vals_type), POINTER :: xb88_params
641 :
642 : CHARACTER(len=*), PARAMETER :: routineN = 'xb88_lsd_eval'
643 :
644 : INTEGER :: handle, i, ispin, npoints
645 : INTEGER, DIMENSION(2, 3) :: bo
646 : REAL(kind=dp) :: epsilon_rho, sx
647 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
648 3092 : POINTER :: dummy, e_0
649 27828 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
650 55656 : e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, e_rho_rho_rho, &
651 27828 : norm_drho, rho, rho_1_3
652 : TYPE(xc_derivative_type), POINTER :: deriv
653 :
654 3092 : CALL timeset(routineN, handle)
655 :
656 3092 : CALL cite_reference(Becke1988)
657 :
658 3092 : NULLIFY (deriv)
659 9276 : DO i = 1, 2
660 9276 : NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
661 : END DO
662 :
663 3092 : CALL section_vals_val_get(xb88_params, "scale_x", r_val=sx)
664 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
665 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
666 : rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
667 : norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
668 3092 : local_bounds=bo)
669 3092 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
670 :
671 3092 : dummy => rho(1)%array
672 :
673 3092 : e_0 => dummy
674 9276 : DO i = 1, 2
675 6184 : e_rho(i)%array => dummy
676 6184 : e_ndrho(i)%array => dummy
677 6184 : e_rho_rho(i)%array => dummy
678 6184 : e_ndrho_rho(i)%array => dummy
679 6184 : e_ndrho_ndrho(i)%array => dummy
680 6184 : e_rho_rho_rho(i)%array => dummy
681 6184 : e_ndrho_rho_rho(i)%array => dummy
682 6184 : e_ndrho_ndrho_rho(i)%array => dummy
683 9276 : e_ndrho_ndrho_ndrho(i)%array => dummy
684 : END DO
685 :
686 3092 : IF (grad_deriv >= 0) THEN
687 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
688 3092 : allocate_deriv=.TRUE.)
689 3092 : CALL xc_derivative_get(deriv, deriv_data=e_0)
690 : END IF
691 3092 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
692 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
693 2996 : allocate_deriv=.TRUE.)
694 2996 : CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
695 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
696 2996 : allocate_deriv=.TRUE.)
697 2996 : CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
698 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
699 2996 : allocate_deriv=.TRUE.)
700 2996 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
701 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
702 2996 : allocate_deriv=.TRUE.)
703 2996 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
704 : END IF
705 3092 : IF (grad_deriv >= 2 .OR. grad_deriv == -2) THEN
706 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
707 4 : allocate_deriv=.TRUE.)
708 4 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(1)%array)
709 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
710 4 : allocate_deriv=.TRUE.)
711 4 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(2)%array)
712 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
713 4 : allocate_deriv=.TRUE.)
714 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(1)%array)
715 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
716 4 : allocate_deriv=.TRUE.)
717 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(2)%array)
718 : deriv => xc_dset_get_derivative(deriv_set, &
719 4 : [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
720 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(1)%array)
721 : deriv => xc_dset_get_derivative(deriv_set, &
722 4 : [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
723 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(2)%array)
724 : END IF
725 3092 : IF (grad_deriv >= 3 .OR. grad_deriv == -3) THEN
726 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
727 0 : allocate_deriv=.TRUE.)
728 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(1)%array)
729 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
730 0 : allocate_deriv=.TRUE.)
731 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(2)%array)
732 : deriv => xc_dset_get_derivative(deriv_set, &
733 0 : [deriv_norm_drhoa, deriv_rhoa, deriv_rhoa], allocate_deriv=.TRUE.)
734 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(1)%array)
735 : deriv => xc_dset_get_derivative(deriv_set, &
736 0 : [deriv_norm_drhob, deriv_rhob, deriv_rhob], allocate_deriv=.TRUE.)
737 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(2)%array)
738 : deriv => xc_dset_get_derivative(deriv_set, &
739 0 : [deriv_norm_drhoa, deriv_norm_drhoa, deriv_rhoa], allocate_deriv=.TRUE.)
740 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(1)%array)
741 : deriv => xc_dset_get_derivative(deriv_set, &
742 0 : [deriv_norm_drhob, deriv_norm_drhob, deriv_rhob], allocate_deriv=.TRUE.)
743 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(2)%array)
744 : deriv => xc_dset_get_derivative(deriv_set, &
745 0 : [deriv_norm_drhoa, deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
746 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(1)%array)
747 : deriv => xc_dset_get_derivative(deriv_set, &
748 0 : [deriv_norm_drhob, deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
749 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(2)%array)
750 : END IF
751 3092 : IF (grad_deriv > 3 .OR. grad_deriv < -3) THEN
752 0 : CPABORT("derivatives bigger than 3 not implemented")
753 : END IF
754 :
755 9276 : DO ispin = 1, 2
756 :
757 : !$OMP PARALLEL DEFAULT(NONE) &
758 : !$OMP SHARED(rho, ispin, rho_1_3, norm_drho, e_0) &
759 : !$OMP SHARED(e_rho, e_ndrho, e_rho_rho, e_ndrho_rho) &
760 : !$OMP SHARED(e_ndrho_ndrho, e_rho_rho_rho) &
761 : !$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
762 : !$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
763 9276 : !$OMP SHARED(epsilon_rho, sx)
764 :
765 : CALL xb88_lsd_calc( &
766 : rho_spin=rho(ispin)%array, &
767 : rho_1_3_spin=rho_1_3(ispin)%array, &
768 : norm_drho_spin=norm_drho(ispin)%array, &
769 : e_0=e_0, &
770 : e_rho_spin=e_rho(ispin)%array, &
771 : e_ndrho_spin=e_ndrho(ispin)%array, &
772 : e_rho_rho_spin=e_rho_rho(ispin)%array, &
773 : e_ndrho_rho_spin=e_ndrho_rho(ispin)%array, &
774 : e_ndrho_ndrho_spin=e_ndrho_ndrho(ispin)%array, &
775 : e_rho_rho_rho_spin=e_rho_rho_rho(ispin)%array, &
776 : e_ndrho_rho_rho_spin=e_ndrho_rho_rho(ispin)%array, &
777 : e_ndrho_ndrho_rho_spin=e_ndrho_ndrho_rho(ispin)%array, &
778 : e_ndrho_ndrho_ndrho_spin=e_ndrho_ndrho_ndrho(ispin)%array, &
779 : grad_deriv=grad_deriv, npoints=npoints, &
780 : epsilon_rho=epsilon_rho, sx=sx)
781 :
782 : !$OMP END PARALLEL
783 :
784 : END DO
785 :
786 3092 : CALL timestop(handle)
787 :
788 3092 : END SUBROUTINE xb88_lsd_eval
789 :
790 : ! **************************************************************************************************
791 : !> \brief low level calculation of the becke 88 exchange functional for lsd
792 : !> \param rho_spin alpha or beta spin density
793 : !> \param rho_1_3_spin rho_spin**(1./3.)
794 : !> \param norm_drho_spin || grad rho_spin ||
795 : !> \param e_0 adds to it the local value of the functional
796 : !> \param e_rho_spin e_*_spin: derivative of the functional wrt. to the variables
797 : !> named where the * is. Everything wrt. to the spin of the arguments.
798 : !> \param e_ndrho_spin ...
799 : !> \param e_rho_rho_spin ...
800 : !> \param e_ndrho_rho_spin ...
801 : !> \param e_ndrho_ndrho_spin ...
802 : !> \param e_rho_rho_rho_spin ...
803 : !> \param e_ndrho_rho_rho_spin ...
804 : !> \param e_ndrho_ndrho_rho_spin ...
805 : !> \param e_ndrho_ndrho_ndrho_spin ...
806 : !> \param grad_deriv ...
807 : !> \param npoints ...
808 : !> \param epsilon_rho ...
809 : !> \param sx scaling-parameter for exchange
810 : !> \par History
811 : !> 01.2004 created [fawzi]
812 : !> 01.2007 added scaling [Manuel Guidon]
813 : !> \author fawzi
814 : ! **************************************************************************************************
815 6184 : SUBROUTINE xb88_lsd_calc(rho_spin, rho_1_3_spin, norm_drho_spin, e_0, &
816 : e_rho_spin, e_ndrho_spin, e_rho_rho_spin, e_ndrho_rho_spin, &
817 : e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
818 : e_ndrho_ndrho_rho_spin, &
819 : e_ndrho_ndrho_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx)
820 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho_spin, rho_1_3_spin, norm_drho_spin
821 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin, e_rho_rho_spin, &
822 : e_ndrho_rho_spin, e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
823 : e_ndrho_ndrho_rho_spin, e_ndrho_ndrho_ndrho_spin
824 : INTEGER, INTENT(in) :: grad_deriv, npoints
825 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
826 :
827 : INTEGER :: ii
828 : REAL(kind=dp) :: c, epsilon_rho43, my_epsilon_rho, my_rho, t1, t103, t11, t12, t127, t133, &
829 : t138, t14, t151, t17, t18, t2, t20, t22, t23, t27, t28, t3, t30, t35, t36, t4, t42, t43, &
830 : t44, t5, t51, t53, t57, t58, t59, t6, t63, t64, t66, t67, t7, t75, t76, t79, t8, t87, x
831 :
832 6184 : c = 1.5_dp*(0.75_dp/Pi)**(1.0_dp/3.0_dp)
833 6184 : my_epsilon_rho = 0.5_dp*epsilon_rho
834 : epsilon_rho43 = my_epsilon_rho**(4._dp/3._dp)
835 :
836 6376 : SELECT CASE (grad_deriv)
837 : CASE (0)
838 :
839 192 : !$OMP DO
840 :
841 : DO ii = 1, npoints
842 14458000 : my_rho = rho_spin(ii)
843 14458000 : IF (my_rho > my_epsilon_rho) THEN
844 14142740 : t1 = rho_1_3_spin(ii)*my_rho
845 14142740 : x = norm_drho_spin(ii)/t1
846 14142740 : t2 = x**2
847 14142740 : t3 = beta*t2
848 14142740 : t4 = beta*x
849 14142740 : t5 = t2 + 0.1e1_dp
850 14142740 : t6 = SQRT(t5)
851 14142740 : t7 = x + t6
852 14142740 : t8 = LOG(t7)
853 14142740 : t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
854 14142740 : t12 = 0.1e1_dp/t11
855 14142740 : t14 = -c - t3*t12
856 :
857 : e_0(ii) = e_0(ii) &
858 14142740 : + (t1*t14)*sx
859 : END IF
860 : END DO
861 :
862 : !$OMP END DO
863 :
864 : CASE (1)
865 :
866 5984 : !$OMP DO
867 :
868 : DO ii = 1, npoints
869 98584998 : my_rho = rho_spin(ii)
870 98584998 : IF (my_rho > my_epsilon_rho) THEN
871 93964520 : t1 = rho_1_3_spin(ii)*my_rho
872 93964520 : x = norm_drho_spin(ii)/t1
873 93964520 : t2 = x**2
874 93964520 : t3 = beta*t2
875 93964520 : t4 = beta*x
876 93964520 : t5 = t2 + 0.1e1_dp
877 93964520 : t6 = SQRT(t5)
878 93964520 : t7 = x + t6
879 93964520 : t8 = LOG(t7)
880 93964520 : t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
881 93964520 : t12 = 0.1e1_dp/t11
882 93964520 : t14 = -c - t3*t12
883 :
884 : e_0(ii) = e_0(ii) &
885 93964520 : + (t1*t14)*sx
886 :
887 93964520 : t17 = t11**2
888 93964520 : t18 = 0.1e1_dp/t17
889 93964520 : t20 = 0.1e1_dp/t6
890 93964520 : t22 = 0.1e1_dp + t20*x
891 93964520 : t23 = 0.1e1_dp/t7
892 93964520 : t27 = 0.6e1_dp*beta*t8 + 0.6e1_dp*t4*t22*t23
893 93964520 : t28 = t18*t27
894 93964520 : t30 = -0.2e1_dp*t4*t12 + t3*t28
895 :
896 : e_rho_spin(ii) = e_rho_spin(ii) &
897 : - (0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t30*x - &
898 93964520 : 0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t14)*sx
899 : e_ndrho_spin(ii) = e_ndrho_spin(ii) &
900 93964520 : + (t30)*sx
901 : END IF
902 : END DO
903 :
904 : !$OMP END DO
905 :
906 : CASE default
907 :
908 6184 : !$OMP DO
909 :
910 : DO ii = 1, npoints
911 159250 : my_rho = rho_spin(ii)
912 159250 : IF (my_rho > my_epsilon_rho) THEN
913 159250 : t1 = rho_1_3_spin(ii)*my_rho
914 159250 : x = norm_drho_spin(ii)/t1
915 159250 : t2 = x**2
916 159250 : t3 = beta*t2
917 159250 : t4 = beta*x
918 159250 : t5 = t2 + 0.1e1_dp
919 159250 : t6 = SQRT(t5)
920 159250 : t7 = x + t6
921 159250 : t8 = LOG(t7)
922 159250 : t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
923 159250 : t12 = 0.1e1_dp/t11
924 159250 : t14 = -c - t3*t12
925 :
926 159250 : IF (grad_deriv >= 0) THEN
927 : e_0(ii) = e_0(ii) &
928 159250 : + (t1*t14)*sx
929 : END IF
930 :
931 159250 : t17 = t11**2
932 159250 : t18 = 0.1e1_dp/t17
933 159250 : t20 = 0.1e1_dp/t6
934 159250 : t22 = 0.1e1_dp + t20*x
935 159250 : t23 = 0.1e1_dp/t7
936 159250 : t27 = 0.6e1_dp*beta*t8 + 0.6e1_dp*t4*t22*t23
937 159250 : t28 = t18*t27
938 159250 : t30 = -0.2e1_dp*t4*t12 + t3*t28
939 :
940 159250 : IF (grad_deriv == -1 .OR. grad_deriv >= 1) THEN
941 : e_rho_spin(ii) = e_rho_spin(ii) &
942 : - (0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t30*x - &
943 159250 : 0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t14)*sx
944 : e_ndrho_spin(ii) = e_ndrho_spin(ii) &
945 159250 : + (t30)*sx
946 : END IF
947 :
948 159250 : t35 = rho_1_3_spin(ii)**2
949 159250 : t36 = 0.1e1_dp/t35
950 159250 : t42 = 0.1e1_dp/t17/t11
951 159250 : t43 = t27**2
952 159250 : t44 = t42*t43
953 159250 : t51 = 0.1e1_dp/t6/t5
954 159250 : t53 = -t51*t2 + t20
955 159250 : t57 = t22**2
956 159250 : t58 = t7**2
957 159250 : t59 = 0.1e1_dp/t58
958 159250 : t63 = 0.12e2_dp*beta*t22*t23 + 0.6e1_dp*t4*t53*t23 - 0.6e1_dp*t4*t57*t59
959 159250 : t64 = t18*t63
960 159250 : t66 = -0.2e1_dp*beta*t12 + 0.4e1_dp*t4*t28 - 0.2e1_dp*t3*t44 + t3*t64
961 159250 : t67 = t36*t66
962 159250 : t75 = 0.1e1_dp/my_rho
963 159250 : t76 = t75*t66
964 159250 : t79 = 0.1e1_dp/t1
965 :
966 159250 : IF (grad_deriv == -2 .OR. grad_deriv >= 2) THEN
967 : e_rho_rho_spin(ii) = e_rho_rho_spin(ii) &
968 : + (0.16e2_dp/0.9e1_dp*t67*t2 - 0.4e1_dp/0.9e1_dp &
969 159250 : *t36*t30*x + 0.4e1_dp/0.9e1_dp*t36*t14)*sx
970 : e_ndrho_rho_spin(ii) = e_ndrho_rho_spin(ii) &
971 159250 : - (0.4e1_dp/0.3e1_dp*t76*x)*sx
972 : e_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_spin(ii) + &
973 159250 : (t66*t79)*sx
974 : END IF
975 :
976 159250 : t87 = t17**2
977 159250 : t103 = t5**2
978 : t127 = 0.6e1_dp*beta*t18*t27 - 0.12e2_dp*t4*t44 + 0.6e1_dp &
979 : *t4*t64 + 0.6e1_dp*t3/t87*t43*t27 - 0.6e1_dp*t3 &
980 : *t42*t27*t63 + t3*t18*(0.18e2_dp*beta*t53*t23 &
981 : - 0.18e2_dp*beta*t57*t59 + 0.6e1_dp*t4*(0.3e1_dp/ &
982 : t6/t103*t2*x - 0.3e1_dp*t51*x)*t23 - 0.18e2_dp* &
983 159250 : t4*t53*t59*t22 + 0.12e2_dp*t4*t57*t22/t58/t7)
984 : t133 = 0.16e2_dp/0.9e1_dp*t36*t127*t2 + 0.28e2_dp/0.9e1_dp &
985 159250 : *t67*x
986 159250 : t138 = 0.1e1_dp/t35/my_rho
987 159250 : t151 = my_rho**2
988 :
989 159250 : IF (grad_deriv == -3 .OR. grad_deriv >= 3) THEN
990 : e_rho_rho_rho_spin(ii) = e_rho_rho_rho_spin(ii) &
991 : - (0.4e1_dp/0.3e1_dp*t133*x*t75 + 0.32e2_dp/0.27e2_dp &
992 : *t138*t66*t2 - 0.8e1_dp/0.27e2_dp*t138*t30* &
993 0 : x + 0.8e1_dp/0.27e2_dp*t138*t14)*sx
994 : e_ndrho_rho_rho_spin(ii) = e_ndrho_rho_rho_spin(ii) &
995 0 : + (t133*t79)*sx
996 : e_ndrho_ndrho_rho_spin(ii) = e_ndrho_ndrho_rho_spin(ii) &
997 : + ((-0.4e1_dp/0.3e1_dp*t75*t127*x - &
998 0 : 0.4e1_dp/0.3e1_dp*t76)*t79)*sx
999 : e_ndrho_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_ndrho_spin(ii) &
1000 0 : + (t127/t35/t151)*sx
1001 : END IF
1002 : END IF
1003 : END DO
1004 :
1005 : !$OMP END DO
1006 :
1007 : END SELECT
1008 :
1009 6184 : END SUBROUTINE xb88_lsd_calc
1010 :
1011 : END MODULE xc_xbecke88
|
