(git:08825b9)
Loading...
Searching...
No Matches
xc_lyp_adiabatic.F
Go to the documentation of this file.
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 density_scaled Lyp functional when used in adiabatic hybrids.
10!> The energy is given as
11!>
12!> Ec = 2*lambda*Ec(rho/lambda) + lambda^2*d/dlambda(Ec(rho/lambda)),
13!>
14!> where rho/lambda is the scaled density
15!> \par History
16!> 1.2008 created [mguidon]
17!> \author Manuel Guidon
18! **************************************************************************************************
19MODULE xc_lyp_adiabatic
20 USE bibliography, ONLY: lee1988,&
21 cite_reference
22 USE input_section_types, ONLY: section_vals_type,&
23 section_vals_val_get
24 USE kinds, ONLY: dp
25 USE mathconstants, ONLY: pi
26 USE xc_derivative_desc, ONLY: deriv_norm_drho,&
27 deriv_norm_drhoa,&
28 deriv_norm_drhob,&
29 deriv_rho,&
30 deriv_rhoa,&
31 deriv_rhob
32 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
33 xc_dset_get_derivative
34 USE xc_derivative_types, ONLY: xc_derivative_get,&
35 xc_derivative_type
36 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
37 USE xc_rho_set_types, ONLY: xc_rho_set_get,&
38 xc_rho_set_type
39#include "../base/base_uses.f90"
40
41 IMPLICIT NONE
42 PRIVATE
43
44 LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
45 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_lyp_adiabatic'
46 REAL(kind=dp), PARAMETER, PRIVATE :: a = 0.04918_dp, b = 0.132_dp, &
47 c = 0.2533_dp, d = 0.349_dp
48
49 PUBLIC :: lyp_adiabatic_lda_info, lyp_adiabatic_lsd_info, lyp_adiabatic_lda_eval, lyp_adiabatic_lsd_eval
50
51CONTAINS
52
53! **************************************************************************************************
54!> \brief return various information on the functional
55!> \param reference string with the reference of the actual functional
56!> \param shortform string with the shortform of the functional name
57!> \param needs the components needed by this functional are set to
58!> true (does not set the unneeded components to false)
59!> \param max_deriv ...
60!> \par History
61!> 01.2008 created [mguidon]
62!> \author Manuel Guidon
63! **************************************************************************************************
64 SUBROUTINE lyp_adiabatic_lda_info(reference, shortform, needs, max_deriv)
65 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
66 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
67 INTEGER, INTENT(out), OPTIONAL :: max_deriv
68
69 IF (PRESENT(reference)) THEN
70 reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LDA version}"
71 END IF
72 IF (PRESENT(shortform)) THEN
73 shortform = "Lee-Yang-Parr correlation energy functional (LDA)"
74 END IF
75 IF (PRESENT(needs)) THEN
76 needs%rho = .true.
77 needs%rho_1_3 = .true.
78 needs%norm_drho = .true.
79 END IF
80 IF (PRESENT(max_deriv)) max_deriv = 1
81
82 END SUBROUTINE lyp_adiabatic_lda_info
83
84! **************************************************************************************************
85!> \brief return various information on the functional
86!> \param reference string with the reference of the actual functional
87!> \param shortform string with the shortform of the functional name
88!> \param needs the components needed by this functional are set to
89!> true (does not set the unneeded components to false)
90!> \param max_deriv ...
91!> \par History
92!> 01.2008 created [mguidon]
93!> \author Manuel Guidon
94! **************************************************************************************************
95 SUBROUTINE lyp_adiabatic_lsd_info(reference, shortform, needs, max_deriv)
96 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
97 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
98 INTEGER, INTENT(out), OPTIONAL :: max_deriv
99
100 IF (PRESENT(reference)) THEN
101 reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LSD version}"
102 END IF
103 IF (PRESENT(shortform)) THEN
104 shortform = "Lee-Yang-Parr correlation energy functional (LSD)"
105 END IF
106 IF (PRESENT(needs)) THEN
107 needs%rho_spin = .true.
108 needs%norm_drho_spin = .true.
109 needs%norm_drho = .true.
110 END IF
111 IF (PRESENT(max_deriv)) max_deriv = 1
112 END SUBROUTINE lyp_adiabatic_lsd_info
113
114! **************************************************************************************************
115!> \brief ...
116!> \param rho_set ...
117!> \param deriv_set ...
118!> \param grad_deriv ...
119!> \param lyp_adiabatic_params ...
120!> \par History
121!> 01.2008 created [mguidon]
122!> \author Manuel Guidon
123! **************************************************************************************************
124 SUBROUTINE lyp_adiabatic_lda_eval(rho_set, deriv_set, grad_deriv, lyp_adiabatic_params)
125 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
126 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
127 INTEGER, INTENT(in) :: grad_deriv
128 TYPE(section_vals_type), POINTER :: lyp_adiabatic_params
129
130 CHARACTER(len=*), PARAMETER :: routinen = 'lyp_adiabatic_lda_eval'
131
132 INTEGER :: handle, npoints
133 INTEGER, DIMENSION(2, 3) :: bo
134 REAL(kind=dp) :: epsilon_norm_drho, epsilon_rho, lambda
135 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
136 POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
137 rho, rho_1_3
138 TYPE(xc_derivative_type), POINTER :: deriv
139
140 CALL timeset(routinen, handle)
141
142 CALL section_vals_val_get(lyp_adiabatic_params, "LAMBDA", r_val=lambda)
143 CALL cite_reference(lee1988)
144
145 CALL xc_rho_set_get(rho_set, rho_1_3=rho_1_3, rho=rho, &
146 norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
147 drho_cutoff=epsilon_norm_drho)
148 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
149
150 dummy => rho
151
152 e_0 => dummy
153 e_rho => dummy
154 e_ndrho => dummy
155
156 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
157 allocate_deriv=.true.)
158 CALL xc_derivative_get(deriv, deriv_data=e_0)
159 IF (grad_deriv >= 1) THEN
160 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
161 allocate_deriv=.true.)
162 CALL xc_derivative_get(deriv, deriv_data=e_rho)
163 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
164 allocate_deriv=.true.)
165 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
166 END IF
167 IF (grad_deriv > 1) THEN
168 cpabort("derivatives bigger than 1 not implemented")
169 END IF
170
171!$OMP PARALLEL DEFAULT(NONE) &
172!$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
173!$OMP SHARED(grad_deriv, npoints) &
174!$OMP SHARED(epsilon_rho, lambda)
175
176 CALL lyp_adiabatic_lda_calc(rho=rho, norm_drho=norm_drho, &
177 e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
178 grad_deriv=grad_deriv, &
179 npoints=npoints, epsilon_rho=epsilon_rho, lambda=lambda)
180
181!$OMP END PARALLEL
182
183 NULLIFY (dummy)
184
185 CALL timestop(handle)
186 END SUBROUTINE lyp_adiabatic_lda_eval
187
188! **************************************************************************************************
189!> \brief ...
190!> \param rho ...
191!> \param norm_drho ...
192!> \param e_0 ...
193!> \param e_rho ...
194!> \param e_ndrho ...
195!> \param grad_deriv ...
196!> \param npoints ...
197!> \param epsilon_rho ...
198!> \param lambda ...
199!> \par History
200!> 01.2008 created [mguidon]
201!> \author Manuel Guidon
202! **************************************************************************************************
203 SUBROUTINE lyp_adiabatic_lda_calc(rho, norm_drho, &
204 e_0, e_rho, e_ndrho, &
205 grad_deriv, npoints, epsilon_rho, lambda)
206 INTEGER, INTENT(in) :: npoints, grad_deriv
207 REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho, e_rho, e_0
208 REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho
209 REAL(kind=dp), INTENT(in) :: epsilon_rho, lambda
210
211 INTEGER :: ii
212 REAL(kind=dp) :: cf, my_ndrho, my_rho, t10, t107, t11, t117, t12, t122, t125, t13, t14, t15, &
213 t153, t16, t17, t180, t189, t19, t195, t2, t20, t25, t28, t29, t3, t34, t36, t37, t38, &
214 t4, t40, t41, t42, t43, t45, t46, t47, t50, t51, t52, t57, t58, t59, t6, t63, t65, t7, &
215 t71, t77, t78, t87, t9, t94
216
217 cf = 0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp)
218
219!$OMP DO
220
221 DO ii = 1, npoints
222 my_rho = rho(ii)
223 IF (my_rho > epsilon_rho) THEN
224 my_ndrho = norm_drho(ii)
225 t2 = d*lambda
226 t3 = my_rho**(0.1e1_dp/0.3e1_dp)
227 t4 = 0.1e1_dp/t3
228 t6 = 0.10e1_dp + t2*t4
229 t7 = 0.1e1_dp/t6
230 t9 = a*b
231 t10 = t9*my_rho
232 t11 = c*lambda
233 t12 = t11*t4
234 t13 = exp(-t12)
235 t14 = t13*t7
236 t15 = my_ndrho**2
237 t16 = my_rho**2
238 t17 = t3**2
239 t19 = 0.1e1_dp/t17/t16
240 t20 = t15*t19
241 t25 = 0.30e1_dp + 0.70e1_dp*t12 + 0.70e1_dp*t2*t4*t7
242 t28 = cf - 0.1388888889e-1_dp*t20*t25
243 t29 = t14*t28
244 t34 = lambda**2
245 t36 = t6**2
246 t37 = 0.1e1_dp/t36
247 t38 = t37*d
248 t40 = t9*t17
249 t41 = c*t13
250 t42 = t7*t28
251 t43 = t41*t42
252 t45 = t13*t37
253 t46 = t28*d
254 t47 = t45*t46
255 t50 = 0.1e1_dp/t17/my_rho
256 t51 = t9*t50
257 t52 = c*t4
258 t57 = d**2
259 t58 = t57*lambda
260 t59 = 0.1e1_dp/t17
261 t63 = 0.70e1_dp*t52 + 0.70e1_dp*d*t4*t7 - 0.70e1_dp*t58*t59*t37
262 t65 = t14*t15*t63
263
264 e_0(ii) = e_0(ii) + 0.20e1_dp*lambda*(-a*my_rho*t7 - t10*t29) + t34*(a*t17 &
265 *t38 + t40*t43 + t40*t47 + 0.13888888888888888889e-1_dp*t51* &
266 t65)
267
268 IF (grad_deriv >= 1) THEN
269 t71 = a*t4
270 t77 = lambda*t13
271 t78 = t77*t42
272 t87 = t16*my_rho
273 t94 = 0.1e1_dp/t3/my_rho
274 t107 = 0.37037037037037037037e-1_dp*t15/t17/t87*t25 - 0.1388888889e-1_dp &
275 *t20*(-0.2333333333e1_dp*t11*t94 - 0.2333333333e1_dp*t2 &
276 *t94*t7 + 0.23333333333333333333e1_dp*t57*t34*t50*t37)
277 t117 = 0.1e1_dp/t36/t6
278 t122 = t9*t4
279 t125 = c**2
280 t153 = 0.1e1_dp/t87
281 t180 = 0.1e1_dp/t16
282 t189 = 0.2e1_dp/0.3e1_dp*t71*t38 + 0.2e1_dp/0.3e1_dp*a*t59*t117* &
283 t57*lambda + 0.2e1_dp/0.3e1_dp*t122*t43 + t9*t59*t125*t78 &
284 /0.3e1_dp + 0.2e1_dp/0.3e1_dp*t9*t59*c*t45*t46*lambda + t40 &
285 *t41*t7*t107 + 0.2e1_dp/0.3e1_dp*t122*t47 + 0.2e1_dp/0.3e1_dp* &
286 t9*t59*t13*t117*t28*t58 + t40*t45*t107*d - 0.2314814815e-1_dp &
287 *t9*t19*t65 + 0.46296296296296296297e-2_dp*t9*t153 &
288 *c*t77*t7*t15*t63 + 0.46296296296296296297e-2_dp*t9*t153 &
289 *t13*t37*t15*t63*d*lambda + 0.13888888888888888889e-1_dp &
290 *t51*t14*t15*(-0.2333333333e1_dp*c*t94 - 0.2333333333e1_dp* &
291 d*t94*t7 + 0.70000000000000000000e1_dp*t57*t50*t37*lambda &
292 - 0.4666666667e1_dp*t57*d*t34*t180*t117)
293
294 e_rho(ii) = e_rho(ii) + 0.20e1_dp*lambda*(-a*t7 - t71*t38*lambda/0.3e1_dp - t9* &
295 t29 - t9*t52*t78/0.3e1_dp - t9*t4*t13*t37*t28*t2/0.3e1_dp &
296 - t10*t14*t107) + t34*t189
297 t195 = t14*my_ndrho*t25
298
299 e_ndrho(ii) = e_ndrho(ii) + 0.55555555555555555556e-1_dp*lambda*a*b*t50*t195 + t34 &
300 *(-0.2777777778e-1_dp*t9*t180*c*t195 - 0.2777777778e-1_dp*t9 &
301 *t180*t13*t37*my_ndrho*t25*d + 0.27777777777777777778e-1_dp* &
302 t51*t14*my_ndrho*t63)
303
304 END IF
305 END IF
306 END DO
307
308!$OMP END DO
309
310 END SUBROUTINE lyp_adiabatic_lda_calc
311
312! **************************************************************************************************
313!> \brief ...
314!> \param rho_set ...
315!> \param deriv_set ...
316!> \param grad_deriv ...
317!> \param lyp_adiabatic_params ...
318!> \par History
319!> 01.2008 created [fawzi]
320!> \author Manuel Guidon
321! **************************************************************************************************
322 SUBROUTINE lyp_adiabatic_lsd_eval(rho_set, deriv_set, grad_deriv, lyp_adiabatic_params)
323 TYPE(xc_rho_set_type) :: rho_set
324 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
325 INTEGER, INTENT(in) :: grad_deriv
326 TYPE(section_vals_type), POINTER :: lyp_adiabatic_params
327
328 CHARACTER(len=*), PARAMETER :: routinen = 'lyp_adiabatic_lsd_eval'
329
330 INTEGER :: handle, npoints
331 INTEGER, DIMENSION(2, 3) :: bo
332 REAL(kind=dp) :: epsilon_rho, lambda
333 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, &
334 e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndra_rb, e_ndrb, e_ndrb_ndrb, &
335 e_ndrb_ra, e_ndrb_rb, e_ra, e_ra_ra, e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, &
336 norm_drhob, rhoa, rhob
337 TYPE(xc_derivative_type), POINTER :: deriv
338
339 CALL timeset(routinen, handle)
340 NULLIFY (deriv)
341
342 CALL section_vals_val_get(lyp_adiabatic_params, "LAMBDA", r_val=lambda)
343 CALL cite_reference(lee1988)
344
345 CALL xc_rho_set_get(rho_set, &
346 rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
347 norm_drhob=norm_drhob, norm_drho=norm_drho, &
348 rho_cutoff=epsilon_rho, &
349 local_bounds=bo)
350 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
351
352 dummy => rhoa
353 e_0 => dummy
354 e_ra => dummy
355 e_rb => dummy
356 e_ndra_ra => dummy
357 e_ndra_rb => dummy
358 e_ndrb_ra => dummy
359 e_ndrb_rb => dummy
360 e_ndr_ndr => dummy
361 e_ndra_ndra => dummy
362 e_ndrb_ndrb => dummy
363 e_ndr => dummy
364 e_ndra => dummy
365 e_ndrb => dummy
366 e_ra_ra => dummy
367 e_ra_rb => dummy
368 e_rb_rb => dummy
369 e_ndr_ra => dummy
370 e_ndr_rb => dummy
371
372 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
373 allocate_deriv=.true.)
374 CALL xc_derivative_get(deriv, deriv_data=e_0)
375 IF (grad_deriv == 1) THEN
376 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
377 allocate_deriv=.true.)
378 CALL xc_derivative_get(deriv, deriv_data=e_ra)
379 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
380 allocate_deriv=.true.)
381 CALL xc_derivative_get(deriv, deriv_data=e_rb)
382 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
383 allocate_deriv=.true.)
384 CALL xc_derivative_get(deriv, deriv_data=e_ndr)
385 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
386 allocate_deriv=.true.)
387 CALL xc_derivative_get(deriv, deriv_data=e_ndra)
388 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
389 allocate_deriv=.true.)
390 CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
391 END IF
392 IF (grad_deriv > 1) THEN
393 cpabort("derivatives bigger than 1 not implemented")
394 END IF
395
396!$OMP PARALLEL DEFAULT(NONE) &
397!$OMP SHARED(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob) &
398!$OMP SHARED(e_0, e_ra, e_rb, e_ndr, e_ndra, e_ndrb) &
399!$OMP SHARED(grad_deriv, npoints) &
400!$OMP SHARED(epsilon_rho, lambda)
401
402 CALL lyp_adiabatic_lsd_calc( &
403 rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa, &
404 norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb, &
405 e_ndr=e_ndr, &
406 e_ndra=e_ndra, e_ndrb=e_ndrb, &
407 grad_deriv=grad_deriv, npoints=npoints, &
408 epsilon_rho=epsilon_rho, lambda=lambda)
409
410!$OMP END PARALLEL
411
412 CALL timestop(handle)
413 END SUBROUTINE lyp_adiabatic_lsd_eval
414
415! **************************************************************************************************
416!> \brief ...
417!> \param rhoa ...
418!> \param rhob ...
419!> \param norm_drho ...
420!> \param norm_drhoa ...
421!> \param norm_drhob ...
422!> \param e_0 ...
423!> \param e_ra ...
424!> \param e_rb ...
425!> \param e_ndr ...
426!> \param e_ndra ...
427!> \param e_ndrb ...
428!> \param grad_deriv ...
429!> \param npoints ...
430!> \param epsilon_rho ...
431!> \param lambda ...
432!> \par History
433!> 08.2008 created [mguidon]
434!> \author Manuel Guidon
435! **************************************************************************************************
436 SUBROUTINE lyp_adiabatic_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, &
437 e_0, e_ra, e_rb, &
438 e_ndr, &
439 e_ndra, e_ndrb, &
440 grad_deriv, npoints, epsilon_rho, lambda)
441 REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, norm_drhoa, &
442 norm_drhob
443 REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ndr, e_ndra, e_ndrb
444 INTEGER, INTENT(in) :: grad_deriv, npoints
445 REAL(kind=dp), INTENT(in) :: epsilon_rho, lambda
446
447 INTEGER :: ii
448 REAL(kind=dp) :: cf, my_ndrho, my_ndrhoa, my_ndrhob, my_rhoa, my_rhob, t1, t10, t100, t102, &
449 t103, t106, t108, t113, t115, t118, t119, t124, t125, t128, t129, t132, t135, t138, t14, &
450 t140, t141, t143, t145, t146, t15, t151, t153, t157, t16, t162, t165, t169, t17, t171, &
451 t174, t179, t18, t183, t186, t187, t188, t19, t194, t196, t199, t2, t200, t202, t21, &
452 t212, t216, t220, t222, t223, t225, t23, t231, t237, t24, t246, t25, t250, t259, t26, &
453 t266, t27, t270, t273, t276, t28, t280, t285, t288, t294, t3, t30, t300, t307, t31, t316, &
454 t32, t325, t348, t351, t355, t362, t387, t39, t394, t4, t41, t42
455 REAL(kind=dp) :: t421, t46, t47, t48, t49, t5, t51, t55, t58, t6, t62, t63, t65, t67, t7, &
456 t73, t74, t76, t77, t78, t80, t83, t84, t85, t86, t87, t9, t90, t91, t94, t95, t96, t97
457
458 cf = 0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp)
459
460!$OMP DO
461
462 DO ii = 1, npoints
463 my_rhoa = max(rhoa(ii), 0.0_dp)
464 my_rhob = max(rhob(ii), 0.0_dp)
465 IF (my_rhoa + my_rhob > epsilon_rho) THEN
466 my_ndrhoa = norm_drhoa(ii)
467 my_ndrhob = norm_drhob(ii)
468 my_ndrho = norm_drho(ii)
469 t1 = a*my_rhoa
470 t2 = my_rhoa + my_rhob
471 t3 = 0.1e1_dp/t2
472 t4 = my_rhob*t3
473 t5 = d*lambda
474 t6 = t2**(0.1e1_dp/0.3e1_dp)
475 t7 = 0.1e1_dp/t6
476 t9 = 0.10e1_dp + t5*t7
477 t10 = 0.1e1_dp/t9
478 t14 = a*b
479 t15 = c*lambda
480 t16 = t15*t7
481 t17 = exp(-t16)
482 t18 = t14*t17
483 t19 = t2**2
484 t21 = t6**2
485 t23 = 0.1e1_dp/t21/t19/t2
486 t24 = t10*t23
487 t25 = my_rhoa*my_rhob
488 t26 = my_rhoa**2
489 t27 = my_rhoa**(0.1e1_dp/0.3e1_dp)
490 t28 = t27**2
491 t30 = my_rhob**2
492 t31 = my_rhob**(0.1e1_dp/0.3e1_dp)
493 t32 = t31**2
494 t39 = t5*t7*t10
495 t41 = 0.26111111111111111111e1_dp - 0.3888888889e0_dp*t16 - 0.3888888889e0_dp &
496 *t39
497 t42 = my_ndrho**2
498 t46 = 0.25000000000000000000e1_dp - 0.5555555556e-1_dp*t16 - 0.5555555556e-1_dp &
499 *t39
500 t47 = my_ndrhoa**2
501 t48 = my_ndrhob**2
502 t49 = t47 + t48
503 t51 = t16 + t39 - 0.110e2_dp
504 t55 = my_rhoa*t3*t47 + t4*t48
505 t58 = 0.12699208415745595798e2_dp*cf*(t28*t26 + t32*t30) + t41 &
506 *t42 - t46*t49 - 0.1111111111e0_dp*t51*t55
507 t62 = 0.66666666666666666667e0_dp*t19
508 t63 = t62 - t26
509 t65 = t62 - t30
510 t67 = t25*t58 - 0.6666666667e0_dp*t19*t42 + t63*t48 + t65*t47
511 t73 = lambda**2
512 t74 = t1*my_rhob
513 t76 = 0.1e1_dp/t6/t2
514 t77 = t9**2
515 t78 = 0.1e1_dp/t77
516 t80 = t76*t78*d
517 t83 = t14*c
518 t84 = t19**2
519 t85 = 0.1e1_dp/t84
520 t86 = t85*t17
521 t87 = t10*t67
522 t90 = t78*t85
523 t91 = t67*d
524 t94 = t17*t10
525 t95 = t14*t94
526 t96 = t23*my_rhoa
527 t97 = c*t7
528 t100 = d*t7*t10
529 t102 = d**2
530 t103 = t102*lambda
531 t106 = t103/t21*t78
532 t108 = -0.3888888889e0_dp*t97 - 0.3888888889e0_dp*t100 + 0.38888888888888888889e0_dp &
533 *t106
534 t113 = -0.5555555556e-1_dp*t97 - 0.5555555556e-1_dp*t100 + 0.55555555555555555556e-1_dp &
535 *t106
536 t115 = t97 + t100 - t106
537 t118 = t108*t42 - t113*t49 - 0.1111111111e0_dp*t115*t55
538 t119 = my_rhob*t118
539
540 e_0(ii) = e_0(ii) + 0.20e1_dp*lambda*(-0.40e1_dp*t1*t4*t10 - t18*t24*t67) &
541 + t73*(0.40e1_dp*t74*t80 + t83*t86*t87 + t18*t90*t91 - &
542 t95*t96*t119)
543
544 IF (grad_deriv == 1) THEN
545 t124 = a*my_rhob
546 t125 = t3*t10
547 t128 = 0.1e1_dp/t19
548 t129 = my_rhob*t128
549 t132 = 0.40e1_dp*t1*t129*t10
550 t135 = 0.1e1_dp/t6/t19*t78
551 t138 = 0.1333333333e1_dp*t74*t135*t5
552 t140 = t84*t2
553 t141 = 0.1e1_dp/t140
554 t143 = t141*t17*t87
555 t145 = t14*t15*t143/0.3e1_dp
556 t146 = t17*t78
557 t151 = t14*t146*t141*t67*t5/0.3e1_dp
558 t153 = 0.1e1_dp/t21/t84
559 t157 = 0.11e2_dp/0.3e1_dp*t18*t10*t153*t67
560 t162 = t15*t76
561 t165 = t5*t76*t10
562 t169 = 0.1e1_dp/t21/t2
563 t171 = t102*t73*t169*t78
564 t174 = (0.12962962962962962963e0_dp*t162 + 0.12962962962962962963e0_dp &
565 *t165 - 0.1296296296e0_dp*t171)*t42
566 t179 = (0.18518518518518518519e-1_dp*t162 + 0.18518518518518518519e-1_dp &
567 *t165 - 0.1851851852e-1_dp*t171)*t49
568 t183 = 0.1111111111e0_dp*(-t162/0.3e1_dp - t165/0.3e1_dp + t171/0.3e1_dp) &
569 *t55
570 t186 = my_rhoa*t128*t47
571 t187 = t129*t48
572 t188 = t3*t47 - t186 - t187
573 t194 = 0.1333333333e1_dp*t2*t42
574 t196 = 0.13333333333333333333e1_dp*my_rhob
575 t199 = 0.13333333333333333333e1_dp*my_rhoa
576 t200 = t199 + t196
577 t202 = my_rhob*t58 + t25*(0.33864555775321588795e2_dp*cf*t28*my_rhoa &
578 + t174 - t179 - t183 - 0.1111111111e0_dp*t51*t188) - t194 + (-0.6666666667e0_dp &
579 *my_rhoa + t196)*t48 + t200*t47
580 t212 = 0.5333333333e1_dp*t74*t135*d
581 t216 = 0.1e1_dp/t77/t9
582 t220 = 0.26666666666666666667e1_dp*t74/t21/t19*t216*t103
583 t222 = 4*t83*t143
584 t223 = c**2
585 t225 = 0.1e1_dp/t6/t140
586 t231 = t14*t223*t225*lambda*t17*t87/0.3e1_dp
587 t237 = 0.2e1_dp/0.3e1_dp*t14*c*t225*t146*t91*lambda
588 t246 = 0.2e1_dp/0.3e1_dp*t14*t17*t216*t225*t67*t103
589 t250 = 4*t18*t78*t141*t91
590 t259 = t14*t15*t141*t94*t25*t118/0.3e1_dp
591 t266 = t14*t146*t141*t25*t118*d*lambda/0.3e1_dp
592 t270 = 0.11e2_dp/0.3e1_dp*t95*t153*my_rhoa*t119
593 t273 = c*t76
594 t276 = d*t76*t10
595 t280 = t102*t169*t78*lambda
596 t285 = t102*d*t73*t128*t216
597 t288 = (0.12962962962962962963e0_dp*t273 + 0.12962962962962962963e0_dp &
598 *t276 - 0.3888888889e0_dp*t280 + 0.25925925925925925926e0_dp*t285) &
599 *t42
600 t294 = (0.18518518518518518519e-1_dp*t273 + 0.18518518518518518519e-1_dp &
601 *t276 - 0.5555555556e-1_dp*t280 + 0.37037037037037037037e-1_dp*t285) &
602 *t49
603 t300 = 0.1111111111e0_dp*(-t273/0.3e1_dp - t276/0.3e1_dp + t280 - 0.2e1_dp &
604 /0.3e1_dp*t285)*t55
605 t307 = 0.40e1_dp*t124*t80 - t212 + t220 - t222 + t231 + t237 + t83 &
606 *t86*t10*t202 + t246 - t250 + t18*t90*t202*d - t259 - &
607 t266 + t270 - t18*t24*t119 - t95*t96*my_rhob*(t288 - t294 - &
608 t300 - 0.1111111111e0_dp*t115*t188)
609
610 e_ra(ii) = e_ra(ii) + 0.20e1_dp*lambda*(-0.40e1_dp*t124*t125 + t132 - t138 - t145 &
611 - t151 + t157 - t18*t24*t202) + t73*t307
612
613 t316 = -t186 + t3*t48 - t187
614 t325 = my_rhoa*t58 + t25*(0.33864555775321588795e2_dp*cf*t32*my_rhob &
615 + t174 - t179 - t183 - 0.1111111111e0_dp*t51*t316) - t194 + t200 &
616 *t48 + (t199 - 0.6666666667e0_dp*my_rhob)*t47
617 t348 = 0.40e1_dp*t1*t80 - t212 + t220 - t222 + t231 + t237 + t83* &
618 t86*t10*t325 + t246 - t250 + t18*t90*t325*d - t259 - t266 &
619 + t270 - t18*t24*my_rhoa*t118 - t95*t96*my_rhob*(t288 - t294 &
620 - t300 - 0.1111111111e0_dp*t115*t316)
621
622 e_rb(ii) = e_rb(ii) + 0.20e1_dp*lambda*(-0.40e1_dp*t1*t125 + t132 - t138 - t145 - &
623 t151 + t157 - t18*t24*t325) + t73*t348
624
625 t351 = lambda*a*b
626 t355 = t3*my_ndrhoa
627 t362 = t25*(-real(2*t46*my_ndrhoa, dp) - 0.2222222222e0_dp*t51*my_rhoa &
628 *t355) + real(2*t65*my_ndrhoa, dp)
629
630 e_ndra(ii) = e_ndra(ii) - 0.20e1_dp*t351*t94*t23*t362 + t73*(t83*t86*t10* &
631 t362 + t18*t90*t362*d - t95*t96*my_rhob*(-real(2*t113* &
632 my_ndrhoa, dp) - 0.2222222222e0_dp*t115*my_rhoa*t355))
633
634 t387 = t3*my_ndrhob
635 t394 = t25*(-real(2*t46*my_ndrhob, dp) - 0.2222222222e0_dp*t51*my_rhob &
636 *t387) + real(2*t63*my_ndrhob, dp)
637
638 e_ndrb(ii) = e_ndrb(ii) - 0.20e1_dp*t351*t94*t23*t394 + t73*(t83*t86*t10* &
639 t394 + t18*t90*t394*d - t95*t96*my_rhob*(-real(2*t113* &
640 my_ndrhob, dp) - 0.2222222222e0_dp*t115*my_rhob*t387))
641
642 t421 = real(2*t25*t41*my_ndrho, dp) - 0.1333333333e1_dp*real(t19, dp)*real(my_ndrho, dp)
643
644 e_ndr(ii) = e_ndr(ii) - 0.20e1_dp*t351*t94*t23*t421 + t73*(t83*t86*t10*t421 &
645 + t18*t90*t421*d - real(2*t95*t96*my_rhob*t108*my_ndrho, dp))
646
647 END IF
648 END IF
649 END DO
650
651!$OMP END DO
652
653 END SUBROUTINE lyp_adiabatic_lsd_calc
654
655END MODULE xc_lyp_adiabatic
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition bibliography.F:21
bibliography::lee1988
integer, save, public lee1988
Definition bibliography.F:36
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition xc_derivative_desc.F:26
xc_derivative_desc::deriv_norm_drho
integer, parameter, public deriv_norm_drho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
Definition xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:65
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_lyp_adiabatic
Calculates the density_scaled Lyp functional when used in adiabatic hybrids. The energy is given as.
Definition xc_lyp_adiabatic.F:19
xc_lyp_adiabatic::lyp_adiabatic_lda_info
subroutine, public lyp_adiabatic_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_lyp_adiabatic.F:65
xc_lyp_adiabatic::lyp_adiabatic_lsd_info
subroutine, public lyp_adiabatic_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_lyp_adiabatic.F:96
xc_lyp_adiabatic::lyp_adiabatic_lda_eval
subroutine, public lyp_adiabatic_lda_eval(rho_set, deriv_set, grad_deriv, lyp_adiabatic_params)
...
Definition xc_lyp_adiabatic.F:125
xc_lyp_adiabatic::lyp_adiabatic_lsd_eval
subroutine, public lyp_adiabatic_lsd_eval(rho_set, deriv_set, grad_deriv, lyp_adiabatic_params)
...
Definition xc_lyp_adiabatic.F:323
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition xc_rho_set_types.F:281
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78