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 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 10624 : 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 10624 : IF (PRESENT(reference)) THEN
64 115 : reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LDA version}"
65 : END IF
66 10624 : IF (PRESENT(shortform)) THEN
67 115 : shortform = "Becke 1988 Exchange Functional (LDA)"
68 : END IF
69 10624 : IF (PRESENT(needs)) THEN
70 10509 : needs%rho = .TRUE.
71 10509 : needs%rho_1_3 = .TRUE.
72 10509 : needs%norm_drho = .TRUE.
73 : END IF
74 10624 : IF (PRESENT(max_deriv)) max_deriv = 3
75 :
76 10624 : 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 2899 : 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 2899 : IF (PRESENT(reference)) THEN
95 39 : reference = "A. Becke, Phys. Rev. A 38, 3098 (1988) {LSD version}"
96 : END IF
97 2899 : IF (PRESENT(shortform)) THEN
98 39 : shortform = "Becke 1988 Exchange Functional (LSD)"
99 : END IF
100 2899 : IF (PRESENT(needs)) THEN
101 2860 : needs%rho_spin = .TRUE.
102 2860 : needs%rho_spin_1_3 = .TRUE.
103 2860 : needs%norm_drho_spin = .TRUE.
104 : END IF
105 2899 : IF (PRESENT(max_deriv)) max_deriv = 3
106 :
107 2899 : 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 maximum degree of the derivative that should be evaluated
115 : !> \param xb88_params input parameters (scaling)
116 : !> \par History
117 : !> 11.2003 created [fawzi]
118 : !> 01.2007 added scaling [Manuel Guidon]
119 : !> \author fawzi
120 : ! **************************************************************************************************
121 28761 : SUBROUTINE xb88_lda_eval(rho_set, deriv_set, grad_deriv, xb88_params)
122 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
123 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
124 : INTEGER, INTENT(in) :: grad_deriv
125 : TYPE(section_vals_type), POINTER :: xb88_params
126 :
127 : CHARACTER(len=*), PARAMETER :: routineN = 'xb88_lda_eval'
128 :
129 : INTEGER :: handle, npoints
130 : INTEGER, DIMENSION(2, 3) :: bo
131 : REAL(kind=dp) :: epsilon_rho, sx
132 9587 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndrho, &
133 9587 : e_ndrho_ndrho, e_ndrho_ndrho_ndrho, e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, &
134 9587 : e_rho, e_rho_rho, e_rho_rho_rho, norm_drho, rho, rho_1_3
135 : TYPE(xc_derivative_type), POINTER :: deriv
136 :
137 9587 : CALL timeset(routineN, handle)
138 :
139 9587 : CALL section_vals_val_get(xb88_params, "scale_x", r_val=sx)
140 :
141 9587 : CALL cite_reference(Becke1988)
142 :
143 : CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
144 9587 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
145 9587 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
146 :
147 9587 : dummy => rho
148 :
149 9587 : e_0 => dummy
150 9587 : e_rho => dummy
151 9587 : e_ndrho => dummy
152 9587 : e_rho_rho => dummy
153 9587 : e_ndrho_rho => dummy
154 9587 : e_ndrho_ndrho => dummy
155 9587 : e_rho_rho_rho => dummy
156 9587 : e_ndrho_rho_rho => dummy
157 9587 : e_ndrho_ndrho_rho => dummy
158 9587 : e_ndrho_ndrho_ndrho => dummy
159 :
160 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
161 9587 : allocate_deriv=.TRUE.)
162 9587 : CALL xc_derivative_get(deriv, deriv_data=e_0)
163 9587 : IF (grad_deriv >= 1) THEN
164 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
165 9343 : allocate_deriv=.TRUE.)
166 9343 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
167 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
168 9343 : allocate_deriv=.TRUE.)
169 9343 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
170 : END IF
171 9587 : IF (grad_deriv >= 2) THEN
172 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
173 372 : allocate_deriv=.TRUE.)
174 372 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
175 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho, deriv_rho], &
176 372 : allocate_deriv=.TRUE.)
177 372 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho)
178 : deriv => xc_dset_get_derivative(deriv_set, &
179 372 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
180 372 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
181 : END IF
182 9587 : IF (grad_deriv >= 3) THEN
183 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
184 0 : allocate_deriv=.TRUE.)
185 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
186 : deriv => xc_dset_get_derivative(deriv_set, &
187 0 : [deriv_norm_drho, deriv_rho, deriv_rho], allocate_deriv=.TRUE.)
188 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho)
189 : deriv => xc_dset_get_derivative(deriv_set, &
190 0 : [deriv_norm_drho, deriv_norm_drho, deriv_rho], allocate_deriv=.TRUE.)
191 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho)
192 : deriv => xc_dset_get_derivative(deriv_set, &
193 0 : [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
194 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
195 : END IF
196 9587 : IF (grad_deriv > 3) THEN
197 0 : CPABORT("derivatives bigger than 3 not implemented")
198 : END IF
199 :
200 : !$OMP PARALLEL DEFAULT(NONE) &
201 : !$OMP SHARED(rho, rho_1_3, norm_drho, e_0, e_rho) &
202 : !$OMP SHARED(e_ndrho, e_rho_rho, e_ndrho_rho) &
203 : !$OMP SHARED(e_ndrho_ndrho, e_rho_rho_rho) &
204 : !$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
205 : !$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
206 9587 : !$OMP SHARED(epsilon_rho, sx)
207 :
208 : CALL xb88_lda_calc(rho=rho, rho_1_3=rho_1_3, norm_drho=norm_drho, &
209 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_rho_rho=e_rho_rho, &
210 : e_ndrho_rho=e_ndrho_rho, e_ndrho_ndrho=e_ndrho_ndrho, &
211 : e_rho_rho_rho=e_rho_rho_rho, e_ndrho_rho_rho=e_ndrho_rho_rho, &
212 : e_ndrho_ndrho_rho=e_ndrho_ndrho_rho, &
213 : e_ndrho_ndrho_ndrho=e_ndrho_ndrho_ndrho, grad_deriv=grad_deriv, &
214 : npoints=npoints, epsilon_rho=epsilon_rho, sx=sx)
215 :
216 : !$OMP END PARALLEL
217 :
218 9587 : CALL timestop(handle)
219 9587 : END SUBROUTINE xb88_lda_eval
220 :
221 : ! **************************************************************************************************
222 : !> \brief evaluates the becke 88 exchange functional for lda
223 : !> \param rho the density where you want to evaluate the functional
224 : !> \param rho_1_3 ...
225 : !> \param norm_drho ...
226 : !> \param e_0 ...
227 : !> \param e_rho ...
228 : !> \param e_ndrho ...
229 : !> \param e_rho_rho ...
230 : !> \param e_ndrho_rho ...
231 : !> \param e_ndrho_ndrho ...
232 : !> \param e_rho_rho_rho ...
233 : !> \param e_ndrho_rho_rho ...
234 : !> \param e_ndrho_ndrho_rho ...
235 : !> \param e_ndrho_ndrho_ndrho ...
236 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
237 : !> \param npoints ...
238 : !> \param epsilon_rho ...
239 : !> \param sx scaling-parameter for exchange
240 : !> \par History
241 : !> 11.2003 created [fawzi]
242 : !> 01.2007 added scaling [Manuel Guidon]
243 : !> \author fawzi
244 : ! **************************************************************************************************
245 9587 : SUBROUTINE xb88_lda_calc(rho, rho_1_3, norm_drho, &
246 9587 : e_0, e_rho, e_ndrho, e_rho_rho, e_ndrho_rho, &
247 9587 : e_ndrho_ndrho, e_rho_rho_rho, e_ndrho_rho_rho, e_ndrho_ndrho_rho, &
248 9587 : e_ndrho_ndrho_ndrho, grad_deriv, npoints, epsilon_rho, sx)
249 : INTEGER, INTENT(in) :: npoints, grad_deriv
250 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho_ndrho_ndrho, &
251 : e_ndrho_ndrho_rho, e_ndrho_rho_rho, e_rho_rho_rho, e_ndrho_ndrho, e_ndrho_rho, e_rho_rho, &
252 : e_ndrho, e_rho, e_0
253 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho_1_3, rho
254 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
255 :
256 : INTEGER :: ii
257 : REAL(kind=dp) :: c, epsilon_rho43, my_rho, my_rho_1_3, t1, t10, t100, t104, t11, t12, t126, &
258 : t13, t159, t164, t17, t170, t176, t18, t189, t2, t21, t23, t29, t3, t31, t33, t37, t39, &
259 : t40, t43, t44, t49, t5, t51, t54, t6, t64, t67, t7, t72, t74, t75, t79, t83, t86, t90, &
260 : t91, t98, t99, x
261 :
262 9587 : c = 1.5_dp*(0.75_dp/Pi)**(1.0_dp/3.0_dp)
263 9587 : t1 = 2**(0.1e1_dp/0.3e1_dp)
264 9587 : t2 = t1**2
265 9587 : t5 = beta*t2
266 9587 : t7 = beta*t1
267 9587 : epsilon_rho43 = epsilon_rho**(4._dp/3._dp)
268 :
269 244 : SELECT CASE (grad_deriv)
270 : CASE (0)
271 :
272 8971 : !$OMP DO
273 :
274 : DO ii = 1, npoints
275 :
276 30441888 : my_rho = rho(ii)
277 30441888 : IF (my_rho > epsilon_rho) THEN
278 24636301 : my_rho_1_3 = rho_1_3(ii)
279 24636301 : t3 = my_rho_1_3*my_rho
280 24636301 : x = norm_drho(ii)/MAX(t3, epsilon_rho43)
281 24636301 : t6 = x**2
282 24636301 : t10 = t2*t6 + 0.1e1_dp
283 24636301 : t11 = SQRT(t10)
284 24636301 : t12 = t1*x + t11
285 24636301 : t13 = LOG(t12)
286 24636301 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
287 24636301 : t18 = 0.1e1_dp/t17
288 24636301 : t21 = -c - t5*t6*t18
289 :
290 24636301 : e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
291 : END IF
292 : END DO
293 :
294 : !$OMP END DO
295 :
296 : CASE (1)
297 :
298 372 : !$OMP DO
299 :
300 : DO ii = 1, npoints
301 :
302 455513332 : my_rho = rho(ii)
303 455513332 : IF (my_rho > epsilon_rho) THEN
304 439052555 : my_rho_1_3 = rho_1_3(ii)
305 439052555 : t3 = my_rho_1_3*my_rho
306 439052555 : x = norm_drho(ii)/t3
307 439052555 : t6 = x**2
308 439052555 : t10 = t2*t6 + 0.1e1_dp
309 439052555 : t11 = SQRT(t10)
310 439052555 : t12 = t1*x + t11
311 439052555 : t13 = LOG(t12)
312 439052555 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
313 439052555 : t18 = 0.1e1_dp/t17
314 439052555 : t21 = -c - t5*t6*t18
315 :
316 439052555 : e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
317 :
318 439052555 : t23 = t2*my_rho_1_3
319 439052555 : t29 = 0.1e1_dp/t11
320 439052555 : t31 = t1 + t2*x*t29
321 439052555 : t33 = 0.1e1_dp/t12
322 439052555 : t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
323 439052555 : t39 = t17**2
324 439052555 : t40 = 0.1e1_dp/t39
325 439052555 : t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
326 439052555 : t44 = t43*x
327 :
328 : e_rho(ii) = e_rho(ii) &
329 : - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
330 439052555 : t23*t21)*sx
331 : e_ndrho(ii) = e_ndrho(ii) &
332 439052555 : + (t2*t43/0.2e1_dp)*sx
333 : END IF
334 :
335 : END DO
336 :
337 : !$OMP END DO
338 :
339 : CASE (2)
340 :
341 0 : !$OMP DO
342 :
343 : DO ii = 1, npoints
344 :
345 5833656 : my_rho = rho(ii)
346 5833656 : IF (my_rho > epsilon_rho) THEN
347 5682366 : my_rho_1_3 = rho_1_3(ii)
348 5682366 : t3 = my_rho_1_3*my_rho
349 5682366 : x = norm_drho(ii)/t3
350 5682366 : t6 = x**2
351 5682366 : t10 = t2*t6 + 0.1e1_dp
352 5682366 : t11 = SQRT(t10)
353 5682366 : t12 = t1*x + t11
354 5682366 : t13 = LOG(t12)
355 5682366 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
356 5682366 : t18 = 0.1e1_dp/t17
357 5682366 : t21 = -c - t5*t6*t18
358 :
359 5682366 : e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
360 :
361 5682366 : t23 = t2*my_rho_1_3
362 5682366 : t29 = 0.1e1_dp/t11
363 5682366 : t31 = t1 + t2*x*t29
364 5682366 : t33 = 0.1e1_dp/t12
365 5682366 : t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
366 5682366 : t39 = t17**2
367 5682366 : t40 = 0.1e1_dp/t39
368 5682366 : t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
369 5682366 : t44 = t43*x
370 :
371 : e_rho(ii) = e_rho(ii) &
372 : - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
373 5682366 : t23*t21)*sx
374 : e_ndrho(ii) = e_ndrho(ii) &
375 5682366 : + (t2*t43/0.2e1_dp)*sx
376 :
377 5682366 : t49 = my_rho_1_3**2
378 5682366 : t51 = t2/t49
379 5682366 : t54 = x*t40
380 5682366 : t64 = 0.1e1_dp/t11/t10
381 5682366 : t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
382 5682366 : t72 = t31**2
383 5682366 : t74 = t12**2
384 5682366 : t75 = 0.1e1_dp/t74
385 : t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
386 5682366 : - 0.6e1_dp*t7*x*t72*t75
387 5682366 : t83 = t37**2
388 5682366 : t86 = 0.1e1_dp/t39/t17
389 : t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
390 5682366 : t79*t40 - 0.2e1_dp*t5*t6*t83*t86
391 5682366 : t91 = t90*t6
392 5682366 : t98 = 0.1e1_dp/my_rho
393 5682366 : t99 = t2*t98
394 5682366 : t100 = t90*x
395 5682366 : t104 = 0.1e1_dp/t3
396 :
397 : e_rho_rho(ii) = e_rho_rho(ii) &
398 : + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
399 5682366 : t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
400 : e_ndrho_rho(ii) = e_ndrho_rho(ii) &
401 5682366 : - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
402 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
403 5682366 : + (t2*t90*t104/0.2e1_dp)*sx
404 : END IF
405 : END DO
406 :
407 : !$OMP END DO
408 :
409 : CASE default
410 :
411 9587 : !$OMP DO
412 :
413 : DO ii = 1, npoints
414 :
415 0 : my_rho = rho(ii)
416 0 : IF (my_rho > epsilon_rho) THEN
417 0 : my_rho_1_3 = rho_1_3(ii)
418 0 : t3 = my_rho_1_3*my_rho
419 0 : x = norm_drho(ii)/t3
420 0 : t6 = x**2
421 0 : t10 = t2*t6 + 0.1e1_dp
422 0 : t11 = SQRT(t10)
423 0 : t12 = t1*x + t11
424 0 : t13 = LOG(t12)
425 0 : t17 = 0.1e1_dp + 0.6e1_dp*t7*x*t13
426 0 : t18 = 0.1e1_dp/t17
427 0 : t21 = -c - t5*t6*t18
428 :
429 0 : e_0(ii) = e_0(ii) + (t2*t3*t21/0.2e1_dp)*sx
430 :
431 0 : IF (grad_deriv >= 1) THEN
432 0 : t23 = t2*my_rho_1_3
433 0 : t29 = 0.1e1_dp/t11
434 0 : t31 = t1 + t2*x*t29
435 0 : t33 = 0.1e1_dp/t12
436 0 : t37 = 0.6e1_dp*t7*t13 + 0.6e1_dp*t7*x*t31*t33
437 0 : t39 = t17**2
438 0 : t40 = 0.1e1_dp/t39
439 0 : t43 = -0.2e1_dp*t5*x*t18 + t5*t6*t37*t40
440 0 : t44 = t43*x
441 :
442 : e_rho(ii) = e_rho(ii) &
443 : - (0.2e1_dp/0.3e1_dp*t23*t44 - 0.2e1_dp/0.3e1_dp* &
444 0 : t23*t21)*sx
445 : e_ndrho(ii) = e_ndrho(ii) &
446 0 : + (t2*t43/0.2e1_dp)*sx
447 : END IF
448 :
449 0 : IF (grad_deriv >= 2) THEN
450 0 : t49 = my_rho_1_3**2
451 0 : t51 = t2/t49
452 0 : t54 = x*t40
453 0 : t64 = 0.1e1_dp/t11/t10
454 0 : t67 = t2*t29 - 0.2e1_dp*t1*t6*t64
455 0 : t72 = t31**2
456 0 : t74 = t12**2
457 0 : t75 = 0.1e1_dp/t74
458 : t79 = 0.12e2_dp*t7*t31*t33 + 0.6e1_dp*t7*x*t67*t33 &
459 0 : - 0.6e1_dp*t7*x*t72*t75
460 0 : t83 = t37**2
461 0 : t86 = 0.1e1_dp/t39/t17
462 : t90 = -0.2e1_dp*t5*t18 + 0.4e1_dp*t5*t54*t37 + t5*t6* &
463 0 : t79*t40 - 0.2e1_dp*t5*t6*t83*t86
464 0 : t91 = t90*t6
465 0 : t98 = 0.1e1_dp/my_rho
466 0 : t99 = t2*t98
467 0 : t100 = t90*x
468 0 : t104 = 0.1e1_dp/t3
469 :
470 : e_rho_rho(ii) = e_rho_rho(ii) &
471 : + (0.8e1_dp/0.9e1_dp*t51*t91 - 0.2e1_dp/0.9e1_dp* &
472 0 : t51*t44 + 0.2e1_dp/0.9e1_dp*t51*t21)*sx
473 : e_ndrho_rho(ii) = e_ndrho_rho(ii) &
474 0 : - (0.2e1_dp/0.3e1_dp*t99*t100)*sx
475 : e_ndrho_ndrho(ii) = e_ndrho_ndrho(ii) &
476 0 : + (t2*t90*t104/0.2e1_dp)*sx
477 : END IF
478 :
479 0 : IF (grad_deriv >= 3) THEN
480 0 : t126 = t10**2
481 0 : t159 = t39**2
482 : t164 = 0.6e1_dp*t5*t40*t37 - 0.12e2_dp*t5*x*t86*t83 &
483 : + 0.6e1_dp*t5*t54*t79 + t5*t6*(0.18e2_dp*t7*t67 &
484 : *t33 - 0.18e2_dp*t7*t72*t75 + 0.6e1_dp*t7*x* &
485 : (-0.6e1_dp*t1*t64*x + 0.12e2_dp*t6*x/t11/t126) &
486 : *t33 - 0.18e2_dp*t7*x*t67*t75*t31 + 0.12e2_dp*t7 &
487 : *x*t72*t31/t74/t12)*t40 - 0.6e1_dp*t5*t6*t79 &
488 0 : *t86*t37 + 0.6e1_dp*t5*t6*t83*t37/t159
489 : t170 = 0.8e1_dp/0.9e1_dp*t51*t164*t6 + 0.14e2_dp/0.9e1_dp &
490 0 : *t51*t100
491 0 : t176 = t2/t49/my_rho
492 0 : t189 = my_rho**2
493 :
494 : e_rho_rho_rho(ii) = e_rho_rho_rho(ii) &
495 : - (0.4e1_dp/0.3e1_dp*t170*x*t98 + 0.16e2_dp/0.27e2_dp &
496 : *t176*t91 - 0.4e1_dp/0.27e2_dp*t176*t44 + &
497 0 : 0.4e1_dp/0.27e2_dp*t176*t21)*sx
498 : e_ndrho_rho_rho(ii) = e_ndrho_rho_rho(ii) &
499 0 : + (t170*t104)*sx
500 : e_ndrho_ndrho_rho(ii) = e_ndrho_ndrho_rho(ii) &
501 : + ((-0.2e1_dp/0.3e1_dp*t99*t164*x - &
502 0 : 0.2e1_dp/0.3e1_dp*t99*t90)*t104)*sx
503 : e_ndrho_ndrho_ndrho(ii) = e_ndrho_ndrho_ndrho(ii) &
504 0 : + (t2*t164/t49/t189/0.2e1_dp)*sx
505 : END IF
506 : END IF
507 : END DO
508 :
509 : !$OMP END DO
510 :
511 : END SELECT
512 :
513 9587 : END SUBROUTINE xb88_lda_calc
514 :
515 : ! **************************************************************************************************
516 : !> \brief evaluates the becke 88 exchange functional for lsd
517 : !> \param rho_set ...
518 : !> \param deriv_set ...
519 : !> \param grad_deriv ...
520 : !> \param xb88_params ...
521 : !> \par History
522 : !> 11.2003 created [fawzi]
523 : !> 01.2007 added scaling [Manuel Guidon]
524 : !> \author fawzi
525 : ! **************************************************************************************************
526 12376 : SUBROUTINE xb88_lsd_eval(rho_set, deriv_set, grad_deriv, xb88_params)
527 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
528 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
529 : INTEGER, INTENT(in) :: grad_deriv
530 : TYPE(section_vals_type), POINTER :: xb88_params
531 :
532 : CHARACTER(len=*), PARAMETER :: routineN = 'xb88_lsd_eval'
533 :
534 : INTEGER :: handle, i, ispin, npoints
535 : INTEGER, DIMENSION(2, 3) :: bo
536 : REAL(kind=dp) :: epsilon_rho, sx
537 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
538 3094 : POINTER :: dummy, e_0
539 27846 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
540 55692 : e_ndrho_ndrho_rho, e_ndrho_rho, e_ndrho_rho_rho, e_rho, e_rho_rho, e_rho_rho_rho, &
541 27846 : norm_drho, rho, rho_1_3
542 : TYPE(xc_derivative_type), POINTER :: deriv
543 :
544 3094 : CALL timeset(routineN, handle)
545 :
546 3094 : CALL cite_reference(Becke1988)
547 :
548 3094 : NULLIFY (deriv)
549 9282 : DO i = 1, 2
550 9282 : NULLIFY (norm_drho(i)%array, rho(i)%array, rho_1_3(i)%array)
551 : END DO
552 :
553 3094 : CALL section_vals_val_get(xb88_params, "scale_x", r_val=sx)
554 : CALL xc_rho_set_get(rho_set, rhoa_1_3=rho_1_3(1)%array, &
555 : rhob_1_3=rho_1_3(2)%array, rhoa=rho(1)%array, &
556 : rhob=rho(2)%array, norm_drhoa=norm_drho(1)%array, &
557 : norm_drhob=norm_drho(2)%array, rho_cutoff=epsilon_rho, &
558 3094 : local_bounds=bo)
559 3094 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
560 :
561 3094 : dummy => rho(1)%array
562 :
563 3094 : e_0 => dummy
564 9282 : DO i = 1, 2
565 6188 : e_rho(i)%array => dummy
566 6188 : e_ndrho(i)%array => dummy
567 6188 : e_rho_rho(i)%array => dummy
568 6188 : e_ndrho_rho(i)%array => dummy
569 6188 : e_ndrho_ndrho(i)%array => dummy
570 6188 : e_rho_rho_rho(i)%array => dummy
571 6188 : e_ndrho_rho_rho(i)%array => dummy
572 6188 : e_ndrho_ndrho_rho(i)%array => dummy
573 9282 : e_ndrho_ndrho_ndrho(i)%array => dummy
574 : END DO
575 :
576 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
577 3094 : allocate_deriv=.TRUE.)
578 3094 : CALL xc_derivative_get(deriv, deriv_data=e_0)
579 3094 : IF (grad_deriv >= 1) THEN
580 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
581 2980 : allocate_deriv=.TRUE.)
582 2980 : CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
583 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
584 2980 : allocate_deriv=.TRUE.)
585 2980 : CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
586 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
587 2980 : allocate_deriv=.TRUE.)
588 2980 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
589 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
590 2980 : allocate_deriv=.TRUE.)
591 2980 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
592 : END IF
593 3094 : IF (grad_deriv >= 2) THEN
594 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa], &
595 4 : allocate_deriv=.TRUE.)
596 4 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(1)%array)
597 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob], &
598 4 : allocate_deriv=.TRUE.)
599 4 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho(2)%array)
600 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa, deriv_rhoa], &
601 4 : allocate_deriv=.TRUE.)
602 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(1)%array)
603 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob, deriv_rhob], &
604 4 : allocate_deriv=.TRUE.)
605 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho(2)%array)
606 : deriv => xc_dset_get_derivative(deriv_set, &
607 4 : [deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
608 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(1)%array)
609 : deriv => xc_dset_get_derivative(deriv_set, &
610 4 : [deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
611 4 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho(2)%array)
612 : END IF
613 3094 : IF (grad_deriv >= 3) THEN
614 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa, deriv_rhoa, deriv_rhoa], &
615 0 : allocate_deriv=.TRUE.)
616 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(1)%array)
617 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob, deriv_rhob, deriv_rhob], &
618 0 : allocate_deriv=.TRUE.)
619 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho(2)%array)
620 : deriv => xc_dset_get_derivative(deriv_set, &
621 0 : [deriv_norm_drhoa, deriv_rhoa, deriv_rhoa], allocate_deriv=.TRUE.)
622 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(1)%array)
623 : deriv => xc_dset_get_derivative(deriv_set, &
624 0 : [deriv_norm_drhob, deriv_rhob, deriv_rhob], allocate_deriv=.TRUE.)
625 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_rho_rho(2)%array)
626 : deriv => xc_dset_get_derivative(deriv_set, &
627 0 : [deriv_norm_drhoa, deriv_norm_drhoa, deriv_rhoa], allocate_deriv=.TRUE.)
628 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(1)%array)
629 : deriv => xc_dset_get_derivative(deriv_set, &
630 0 : [deriv_norm_drhob, deriv_norm_drhob, deriv_rhob], allocate_deriv=.TRUE.)
631 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_rho(2)%array)
632 : deriv => xc_dset_get_derivative(deriv_set, &
633 0 : [deriv_norm_drhoa, deriv_norm_drhoa, deriv_norm_drhoa], allocate_deriv=.TRUE.)
634 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(1)%array)
635 : deriv => xc_dset_get_derivative(deriv_set, &
636 0 : [deriv_norm_drhob, deriv_norm_drhob, deriv_norm_drhob], allocate_deriv=.TRUE.)
637 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho(2)%array)
638 : END IF
639 3094 : IF (grad_deriv > 3) THEN
640 0 : CPABORT("derivatives bigger than 3 not implemented")
641 : END IF
642 :
643 9282 : DO ispin = 1, 2
644 :
645 : !$OMP PARALLEL DEFAULT(NONE) &
646 : !$OMP SHARED(rho, ispin, rho_1_3, norm_drho, e_0) &
647 : !$OMP SHARED(e_rho, e_ndrho, e_rho_rho, e_ndrho_rho) &
648 : !$OMP SHARED(e_ndrho_ndrho, e_rho_rho_rho) &
649 : !$OMP SHARED(e_ndrho_rho_rho, e_ndrho_ndrho_rho) &
650 : !$OMP SHARED(e_ndrho_ndrho_ndrho, grad_deriv, npoints) &
651 9282 : !$OMP SHARED(epsilon_rho, sx)
652 :
653 : CALL xb88_lsd_calc( &
654 : rho_spin=rho(ispin)%array, &
655 : rho_1_3_spin=rho_1_3(ispin)%array, &
656 : norm_drho_spin=norm_drho(ispin)%array, &
657 : e_0=e_0, &
658 : e_rho_spin=e_rho(ispin)%array, &
659 : e_ndrho_spin=e_ndrho(ispin)%array, &
660 : e_rho_rho_spin=e_rho_rho(ispin)%array, &
661 : e_ndrho_rho_spin=e_ndrho_rho(ispin)%array, &
662 : e_ndrho_ndrho_spin=e_ndrho_ndrho(ispin)%array, &
663 : e_rho_rho_rho_spin=e_rho_rho_rho(ispin)%array, &
664 : e_ndrho_rho_rho_spin=e_ndrho_rho_rho(ispin)%array, &
665 : e_ndrho_ndrho_rho_spin=e_ndrho_ndrho_rho(ispin)%array, &
666 : e_ndrho_ndrho_ndrho_spin=e_ndrho_ndrho_ndrho(ispin)%array, &
667 : grad_deriv=grad_deriv, npoints=npoints, &
668 : epsilon_rho=epsilon_rho, sx=sx)
669 :
670 : !$OMP END PARALLEL
671 :
672 : END DO
673 :
674 3094 : CALL timestop(handle)
675 :
676 3094 : END SUBROUTINE xb88_lsd_eval
677 :
678 : ! **************************************************************************************************
679 : !> \brief low level calculation of the becke 88 exchange functional for lsd
680 : !> \param rho_spin alpha or beta spin density
681 : !> \param rho_1_3_spin rho_spin**(1./3.)
682 : !> \param norm_drho_spin || grad rho_spin ||
683 : !> \param e_0 adds to it the local value of the functional
684 : !> \param e_rho_spin e_*_spin: derivative of the functional wrt. to the variables
685 : !> named where the * is. Everything wrt. to the spin of the arguments.
686 : !> \param e_ndrho_spin ...
687 : !> \param e_rho_rho_spin ...
688 : !> \param e_ndrho_rho_spin ...
689 : !> \param e_ndrho_ndrho_spin ...
690 : !> \param e_rho_rho_rho_spin ...
691 : !> \param e_ndrho_rho_rho_spin ...
692 : !> \param e_ndrho_ndrho_rho_spin ...
693 : !> \param e_ndrho_ndrho_ndrho_spin ...
694 : !> \param grad_deriv ...
695 : !> \param npoints ...
696 : !> \param epsilon_rho ...
697 : !> \param sx scaling-parameter for exchange
698 : !> \par History
699 : !> 01.2004 created [fawzi]
700 : !> 01.2007 added scaling [Manuel Guidon]
701 : !> \author fawzi
702 : ! **************************************************************************************************
703 6188 : SUBROUTINE xb88_lsd_calc(rho_spin, rho_1_3_spin, norm_drho_spin, e_0, &
704 : e_rho_spin, e_ndrho_spin, e_rho_rho_spin, e_ndrho_rho_spin, &
705 : e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
706 : e_ndrho_ndrho_rho_spin, &
707 : e_ndrho_ndrho_ndrho_spin, grad_deriv, npoints, epsilon_rho, sx)
708 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho_spin, rho_1_3_spin, norm_drho_spin
709 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho_spin, e_ndrho_spin, e_rho_rho_spin, &
710 : e_ndrho_rho_spin, e_ndrho_ndrho_spin, e_rho_rho_rho_spin, e_ndrho_rho_rho_spin, &
711 : e_ndrho_ndrho_rho_spin, e_ndrho_ndrho_ndrho_spin
712 : INTEGER, INTENT(in) :: grad_deriv, npoints
713 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx
714 :
715 : INTEGER :: ii
716 : REAL(kind=dp) :: c, epsilon_rho43, my_epsilon_rho, my_rho, t1, t103, t11, t12, t127, t133, &
717 : t138, t14, t151, t17, t18, t2, t20, t22, t23, t27, t28, t3, t30, t35, t36, t4, t42, t43, &
718 : t44, t5, t51, t53, t57, t58, t59, t6, t63, t64, t66, t67, t7, t75, t76, t79, t8, t87, x
719 :
720 6188 : c = 1.5_dp*(0.75_dp/Pi)**(1.0_dp/3.0_dp)
721 6188 : my_epsilon_rho = 0.5_dp*epsilon_rho
722 : epsilon_rho43 = my_epsilon_rho**(4._dp/3._dp)
723 :
724 6416 : SELECT CASE (grad_deriv)
725 : CASE (0)
726 :
727 228 : !$OMP DO
728 :
729 : DO ii = 1, npoints
730 16960400 : my_rho = rho_spin(ii)
731 16960400 : IF (my_rho > my_epsilon_rho) THEN
732 16942527 : t1 = rho_1_3_spin(ii)*my_rho
733 16942527 : x = norm_drho_spin(ii)/t1
734 16942527 : t2 = x**2
735 16942527 : t3 = beta*t2
736 16942527 : t4 = beta*x
737 16942527 : t5 = t2 + 0.1e1_dp
738 16942527 : t6 = SQRT(t5)
739 16942527 : t7 = x + t6
740 16942527 : t8 = LOG(t7)
741 16942527 : t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
742 16942527 : t12 = 0.1e1_dp/t11
743 16942527 : t14 = -c - t3*t12
744 :
745 : e_0(ii) = e_0(ii) &
746 16942527 : + (t1*t14)*sx
747 : END IF
748 : END DO
749 :
750 : !$OMP END DO
751 :
752 : CASE (1)
753 :
754 5952 : !$OMP DO
755 :
756 : DO ii = 1, npoints
757 97576692 : my_rho = rho_spin(ii)
758 97576692 : IF (my_rho > my_epsilon_rho) THEN
759 93557774 : t1 = rho_1_3_spin(ii)*my_rho
760 93557774 : x = norm_drho_spin(ii)/t1
761 93557774 : t2 = x**2
762 93557774 : t3 = beta*t2
763 93557774 : t4 = beta*x
764 93557774 : t5 = t2 + 0.1e1_dp
765 93557774 : t6 = SQRT(t5)
766 93557774 : t7 = x + t6
767 93557774 : t8 = LOG(t7)
768 93557774 : t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
769 93557774 : t12 = 0.1e1_dp/t11
770 93557774 : t14 = -c - t3*t12
771 :
772 : e_0(ii) = e_0(ii) &
773 93557774 : + (t1*t14)*sx
774 :
775 93557774 : t17 = t11**2
776 93557774 : t18 = 0.1e1_dp/t17
777 93557774 : t20 = 0.1e1_dp/t6
778 93557774 : t22 = 0.1e1_dp + t20*x
779 93557774 : t23 = 0.1e1_dp/t7
780 93557774 : t27 = 0.6e1_dp*beta*t8 + 0.6e1_dp*t4*t22*t23
781 93557774 : t28 = t18*t27
782 93557774 : t30 = -0.2e1_dp*t4*t12 + t3*t28
783 :
784 : e_rho_spin(ii) = e_rho_spin(ii) &
785 : - (0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t30*x - &
786 93557774 : 0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t14)*sx
787 : e_ndrho_spin(ii) = e_ndrho_spin(ii) &
788 93557774 : + (t30)*sx
789 : END IF
790 : END DO
791 :
792 : !$OMP END DO
793 :
794 : CASE default
795 :
796 6188 : !$OMP DO
797 :
798 : DO ii = 1, npoints
799 159250 : my_rho = rho_spin(ii)
800 159250 : IF (my_rho > my_epsilon_rho) THEN
801 159250 : t1 = rho_1_3_spin(ii)*my_rho
802 159250 : x = norm_drho_spin(ii)/t1
803 159250 : t2 = x**2
804 159250 : t3 = beta*t2
805 159250 : t4 = beta*x
806 159250 : t5 = t2 + 0.1e1_dp
807 159250 : t6 = SQRT(t5)
808 159250 : t7 = x + t6
809 159250 : t8 = LOG(t7)
810 159250 : t11 = 0.1e1_dp + 0.6e1_dp*t4*t8
811 159250 : t12 = 0.1e1_dp/t11
812 159250 : t14 = -c - t3*t12
813 :
814 : e_0(ii) = e_0(ii) &
815 159250 : + (t1*t14)*sx
816 :
817 159250 : IF (grad_deriv >= 1) THEN
818 159250 : t17 = t11**2
819 159250 : t18 = 0.1e1_dp/t17
820 159250 : t20 = 0.1e1_dp/t6
821 159250 : t22 = 0.1e1_dp + t20*x
822 159250 : t23 = 0.1e1_dp/t7
823 159250 : t27 = 0.6e1_dp*beta*t8 + 0.6e1_dp*t4*t22*t23
824 159250 : t28 = t18*t27
825 159250 : t30 = -0.2e1_dp*t4*t12 + t3*t28
826 :
827 : e_rho_spin(ii) = e_rho_spin(ii) &
828 : - (0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t30*x - &
829 159250 : 0.4e1_dp/0.3e1_dp*rho_1_3_spin(ii)*t14)*sx
830 : e_ndrho_spin(ii) = e_ndrho_spin(ii) &
831 159250 : + (t30)*sx
832 : END IF
833 :
834 159250 : IF (grad_deriv >= 2) THEN
835 159250 : t35 = rho_1_3_spin(ii)**2
836 159250 : t36 = 0.1e1_dp/t35
837 159250 : t42 = 0.1e1_dp/t17/t11
838 159250 : t43 = t27**2
839 159250 : t44 = t42*t43
840 159250 : t51 = 0.1e1_dp/t6/t5
841 159250 : t53 = -t51*t2 + t20
842 159250 : t57 = t22**2
843 159250 : t58 = t7**2
844 159250 : t59 = 0.1e1_dp/t58
845 159250 : t63 = 0.12e2_dp*beta*t22*t23 + 0.6e1_dp*t4*t53*t23 - 0.6e1_dp*t4*t57*t59
846 159250 : t64 = t18*t63
847 159250 : t66 = -0.2e1_dp*beta*t12 + 0.4e1_dp*t4*t28 - 0.2e1_dp*t3*t44 + t3*t64
848 159250 : t67 = t36*t66
849 159250 : t75 = 0.1e1_dp/my_rho
850 159250 : t76 = t75*t66
851 159250 : t79 = 0.1e1_dp/t1
852 :
853 : e_rho_rho_spin(ii) = e_rho_rho_spin(ii) &
854 : + (0.16e2_dp/0.9e1_dp*t67*t2 - 0.4e1_dp/0.9e1_dp &
855 159250 : *t36*t30*x + 0.4e1_dp/0.9e1_dp*t36*t14)*sx
856 : e_ndrho_rho_spin(ii) = e_ndrho_rho_spin(ii) &
857 159250 : - (0.4e1_dp/0.3e1_dp*t76*x)*sx
858 : e_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_spin(ii) + &
859 159250 : (t66*t79)*sx
860 : END IF
861 :
862 159250 : IF (grad_deriv >= 3) THEN
863 0 : t87 = t17**2
864 0 : t103 = t5**2
865 : t127 = 0.6e1_dp*beta*t18*t27 - 0.12e2_dp*t4*t44 + 0.6e1_dp &
866 : *t4*t64 + 0.6e1_dp*t3/t87*t43*t27 - 0.6e1_dp*t3 &
867 : *t42*t27*t63 + t3*t18*(0.18e2_dp*beta*t53*t23 &
868 : - 0.18e2_dp*beta*t57*t59 + 0.6e1_dp*t4*(0.3e1_dp/ &
869 : t6/t103*t2*x - 0.3e1_dp*t51*x)*t23 - 0.18e2_dp* &
870 0 : t4*t53*t59*t22 + 0.12e2_dp*t4*t57*t22/t58/t7)
871 : t133 = 0.16e2_dp/0.9e1_dp*t36*t127*t2 + 0.28e2_dp/0.9e1_dp &
872 0 : *t67*x
873 0 : t138 = 0.1e1_dp/t35/my_rho
874 0 : t151 = my_rho**2
875 :
876 : e_rho_rho_rho_spin(ii) = e_rho_rho_rho_spin(ii) &
877 : - (0.4e1_dp/0.3e1_dp*t133*x*t75 + 0.32e2_dp/0.27e2_dp &
878 : *t138*t66*t2 - 0.8e1_dp/0.27e2_dp*t138*t30* &
879 0 : x + 0.8e1_dp/0.27e2_dp*t138*t14)*sx
880 : e_ndrho_rho_rho_spin(ii) = e_ndrho_rho_rho_spin(ii) &
881 0 : + (t133*t79)*sx
882 : e_ndrho_ndrho_rho_spin(ii) = e_ndrho_ndrho_rho_spin(ii) &
883 : + ((-0.4e1_dp/0.3e1_dp*t75*t127*x - &
884 0 : 0.4e1_dp/0.3e1_dp*t76)*t79)*sx
885 : e_ndrho_ndrho_ndrho_spin(ii) = e_ndrho_ndrho_ndrho_spin(ii) &
886 0 : + (t127/t35/t151)*sx
887 : END IF
888 : END IF
889 : END DO
890 :
891 : !$OMP END DO
892 :
893 : END SELECT
894 :
895 6188 : END SUBROUTINE xb88_lsd_calc
896 :
897 : END MODULE xc_xbecke88
|
