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 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 : ! **************************************************************************************************
19 : MODULE 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 :
51 : CONTAINS
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 109 : 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 109 : IF (PRESENT(reference)) THEN
70 1 : reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LDA version}"
71 : END IF
72 109 : IF (PRESENT(shortform)) THEN
73 1 : shortform = "Lee-Yang-Parr correlation energy functional (LDA)"
74 : END IF
75 109 : IF (PRESENT(needs)) THEN
76 108 : needs%rho = .TRUE.
77 108 : needs%rho_1_3 = .TRUE.
78 108 : needs%norm_drho = .TRUE.
79 : END IF
80 109 : IF (PRESENT(max_deriv)) max_deriv = 1
81 :
82 109 : 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 125 : 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 125 : IF (PRESENT(reference)) THEN
101 1 : reference = "C. Lee, W. Yang, R.G. Parr, Phys. Rev. B, 37, 785 (1988) {LSD version}"
102 : END IF
103 125 : IF (PRESENT(shortform)) THEN
104 1 : shortform = "Lee-Yang-Parr correlation energy functional (LSD)"
105 : END IF
106 125 : IF (PRESENT(needs)) THEN
107 124 : needs%rho_spin = .TRUE.
108 124 : needs%norm_drho_spin = .TRUE.
109 124 : needs%norm_drho = .TRUE.
110 : END IF
111 125 : IF (PRESENT(max_deriv)) max_deriv = 1
112 125 : 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 312 : 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 104 : POINTER :: dummy, e_0, e_ndrho, e_rho, norm_drho, &
137 104 : rho, rho_1_3
138 : TYPE(xc_derivative_type), POINTER :: deriv
139 :
140 104 : CALL timeset(routineN, handle)
141 :
142 104 : CALL section_vals_val_get(lyp_adiabatic_params, "LAMBDA", r_val=lambda)
143 104 : 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 104 : drho_cutoff=epsilon_norm_drho)
148 104 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
149 :
150 104 : dummy => rho
151 :
152 104 : e_0 => dummy
153 104 : e_rho => dummy
154 104 : e_ndrho => dummy
155 :
156 104 : IF (grad_deriv >= 0) THEN
157 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
158 104 : allocate_deriv=.TRUE.)
159 104 : CALL xc_derivative_get(deriv, deriv_data=e_0)
160 : END IF
161 104 : IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
162 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
163 66 : allocate_deriv=.TRUE.)
164 66 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
165 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
166 66 : allocate_deriv=.TRUE.)
167 66 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
168 : END IF
169 104 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
170 0 : CPABORT("derivatives bigger than 1 not implemented")
171 : END IF
172 :
173 : !$OMP PARALLEL DEFAULT(NONE) &
174 : !$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
175 : !$OMP SHARED(grad_deriv, npoints) &
176 104 : !$OMP SHARED(epsilon_rho, lambda)
177 :
178 : CALL lyp_adiabatic_lda_calc(rho=rho, norm_drho=norm_drho, &
179 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
180 : grad_deriv=grad_deriv, &
181 : npoints=npoints, epsilon_rho=epsilon_rho, lambda=lambda)
182 :
183 : !$OMP END PARALLEL
184 :
185 104 : NULLIFY (dummy)
186 :
187 104 : CALL timestop(handle)
188 104 : END SUBROUTINE lyp_adiabatic_lda_eval
189 :
190 : ! **************************************************************************************************
191 : !> \brief ...
192 : !> \param rho ...
193 : !> \param norm_drho ...
194 : !> \param e_0 ...
195 : !> \param e_rho ...
196 : !> \param e_ndrho ...
197 : !> \param grad_deriv ...
198 : !> \param npoints ...
199 : !> \param epsilon_rho ...
200 : !> \param lambda ...
201 : !> \par History
202 : !> 01.2008 created [mguidon]
203 : !> \author Manuel Guidon
204 : ! **************************************************************************************************
205 104 : SUBROUTINE lyp_adiabatic_lda_calc(rho, norm_drho, &
206 104 : e_0, e_rho, e_ndrho, &
207 : grad_deriv, npoints, epsilon_rho, lambda)
208 : INTEGER, INTENT(in) :: npoints, grad_deriv
209 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: e_ndrho, e_rho, e_0
210 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(in) :: norm_drho, rho
211 : REAL(kind=dp), INTENT(in) :: epsilon_rho, lambda
212 :
213 : INTEGER :: ii
214 : REAL(kind=dp) :: cf, my_ndrho, my_rho, t10, t107, t11, t117, t12, t122, t125, t13, t14, t15, &
215 : t153, t16, t17, t180, t189, t19, t195, t2, t20, t25, t28, t29, t3, t34, t36, t37, t38, &
216 : t4, t40, t41, t42, t43, t45, t46, t47, t50, t51, t52, t57, t58, t59, t6, t63, t65, t7, &
217 : t71, t77, t78, t87, t9, t94
218 :
219 104 : cf = 0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp)
220 :
221 : !$OMP DO
222 :
223 : DO ii = 1, npoints
224 4224064 : my_rho = rho(ii)
225 4224064 : IF (my_rho > epsilon_rho) THEN
226 1289946 : IF (grad_deriv >= 0) THEN
227 1289946 : my_ndrho = norm_drho(ii)
228 1289946 : t2 = d*lambda
229 1289946 : t3 = my_rho**(0.1e1_dp/0.3e1_dp)
230 1289946 : t4 = 0.1e1_dp/t3
231 1289946 : t6 = 0.10e1_dp + t2*t4
232 1289946 : t7 = 0.1e1_dp/t6
233 1289946 : t9 = a*b
234 1289946 : t10 = t9*my_rho
235 1289946 : t11 = c*lambda
236 1289946 : t12 = t11*t4
237 1289946 : t13 = EXP(-t12)
238 1289946 : t14 = t13*t7
239 1289946 : t15 = my_ndrho**2
240 1289946 : t16 = my_rho**2
241 1289946 : t17 = t3**2
242 1289946 : t19 = 0.1e1_dp/t17/t16
243 1289946 : t20 = t15*t19
244 1289946 : t25 = 0.30e1_dp + 0.70e1_dp*t12 + 0.70e1_dp*t2*t4*t7
245 1289946 : t28 = Cf - 0.1388888889e-1_dp*t20*t25
246 1289946 : t29 = t14*t28
247 1289946 : t34 = lambda**2
248 1289946 : t36 = t6**2
249 1289946 : t37 = 0.1e1_dp/t36
250 1289946 : t38 = t37*d
251 1289946 : t40 = t9*t17
252 1289946 : t41 = c*t13
253 1289946 : t42 = t7*t28
254 1289946 : t43 = t41*t42
255 1289946 : t45 = t13*t37
256 1289946 : t46 = t28*d
257 1289946 : t47 = t45*t46
258 1289946 : t50 = 0.1e1_dp/t17/my_rho
259 1289946 : t51 = t9*t50
260 1289946 : t52 = c*t4
261 1289946 : t57 = d**2
262 1289946 : t58 = t57*lambda
263 1289946 : t59 = 0.1e1_dp/t17
264 1289946 : t63 = 0.70e1_dp*t52 + 0.70e1_dp*d*t4*t7 - 0.70e1_dp*t58*t59*t37
265 1289946 : t65 = t14*t15*t63
266 :
267 : e_0(ii) = e_0(ii) + 0.20e1_dp*lambda*(-a*my_rho*t7 - t10*t29) + t34*(a*t17 &
268 : *t38 + t40*t43 + t40*t47 + 0.13888888888888888889e-1_dp*t51* &
269 1289946 : t65)
270 :
271 : END IF
272 1289946 : IF (grad_deriv >= 1) THEN
273 419495 : t71 = a*t4
274 419495 : t77 = lambda*t13
275 419495 : t78 = t77*t42
276 419495 : t87 = t16*my_rho
277 419495 : t94 = 0.1e1_dp/t3/my_rho
278 : t107 = 0.37037037037037037037e-1_dp*t15/t17/t87*t25 - 0.1388888889e-1_dp &
279 : *t20*(-0.2333333333e1_dp*t11*t94 - 0.2333333333e1_dp*t2 &
280 419495 : *t94*t7 + 0.23333333333333333333e1_dp*t57*t34*t50*t37)
281 419495 : t117 = 0.1e1_dp/t36/t6
282 419495 : t122 = t9*t4
283 419495 : t125 = c**2
284 419495 : t153 = 0.1e1_dp/t87
285 419495 : t180 = 0.1e1_dp/t16
286 : t189 = 0.2e1_dp/0.3e1_dp*t71*t38 + 0.2e1_dp/0.3e1_dp*a*t59*t117* &
287 : t57*lambda + 0.2e1_dp/0.3e1_dp*t122*t43 + t9*t59*t125*t78 &
288 : /0.3e1_dp + 0.2e1_dp/0.3e1_dp*t9*t59*c*t45*t46*lambda + t40 &
289 : *t41*t7*t107 + 0.2e1_dp/0.3e1_dp*t122*t47 + 0.2e1_dp/0.3e1_dp* &
290 : t9*t59*t13*t117*t28*t58 + t40*t45*t107*d - 0.2314814815e-1_dp &
291 : *t9*t19*t65 + 0.46296296296296296297e-2_dp*t9*t153 &
292 : *c*t77*t7*t15*t63 + 0.46296296296296296297e-2_dp*t9*t153 &
293 : *t13*t37*t15*t63*d*lambda + 0.13888888888888888889e-1_dp &
294 : *t51*t14*t15*(-0.2333333333e1_dp*c*t94 - 0.2333333333e1_dp* &
295 : d*t94*t7 + 0.70000000000000000000e1_dp*t57*t50*t37*lambda &
296 419495 : - 0.4666666667e1_dp*t57*d*t34*t180*t117)
297 :
298 : e_rho(ii) = e_rho(ii) + 0.20e1_dp*lambda*(-a*t7 - t71*t38*lambda/0.3e1_dp - t9* &
299 : t29 - t9*t52*t78/0.3e1_dp - t9*t4*t13*t37*t28*t2/0.3e1_dp &
300 419495 : - t10*t14*t107) + t34*t189
301 419495 : t195 = t14*my_ndrho*t25
302 :
303 : e_ndrho(ii) = e_ndrho(ii) + 0.55555555555555555556e-1_dp*lambda*a*b*t50*t195 + t34 &
304 : *(-0.2777777778e-1_dp*t9*t180*c*t195 - 0.2777777778e-1_dp*t9 &
305 : *t180*t13*t37*my_ndrho*t25*d + 0.27777777777777777778e-1_dp* &
306 419495 : t51*t14*my_ndrho*t63)
307 :
308 : END IF
309 : END IF
310 : END DO
311 :
312 : !$OMP END DO
313 :
314 104 : END SUBROUTINE lyp_adiabatic_lda_calc
315 :
316 : ! **************************************************************************************************
317 : !> \brief ...
318 : !> \param rho_set ...
319 : !> \param deriv_set ...
320 : !> \param grad_deriv ...
321 : !> \param lyp_adiabatic_params ...
322 : !> \par History
323 : !> 01.2008 created [fawzi]
324 : !> \author Manuel Guidon
325 : ! **************************************************************************************************
326 360 : SUBROUTINE lyp_adiabatic_lsd_eval(rho_set, deriv_set, grad_deriv, lyp_adiabatic_params)
327 : TYPE(xc_rho_set_type) :: rho_set
328 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
329 : INTEGER, INTENT(in) :: grad_deriv
330 : TYPE(section_vals_type), POINTER :: lyp_adiabatic_params
331 :
332 : CHARACTER(len=*), PARAMETER :: routineN = 'lyp_adiabatic_lsd_eval'
333 :
334 : INTEGER :: handle, npoints
335 : INTEGER, DIMENSION(2, 3) :: bo
336 : REAL(kind=dp) :: epsilon_rho, lambda
337 120 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: dummy, e_0, e_ndr, e_ndr_ndr, &
338 120 : e_ndr_ra, e_ndr_rb, e_ndra, e_ndra_ndra, e_ndra_ra, e_ndra_rb, e_ndrb, e_ndrb_ndrb, &
339 120 : e_ndrb_ra, e_ndrb_rb, e_ra, e_ra_ra, e_ra_rb, e_rb, e_rb_rb, norm_drho, norm_drhoa, &
340 120 : norm_drhob, rhoa, rhob
341 : TYPE(xc_derivative_type), POINTER :: deriv
342 :
343 120 : CALL timeset(routineN, handle)
344 120 : NULLIFY (deriv)
345 :
346 120 : CALL section_vals_val_get(lyp_adiabatic_params, "LAMBDA", r_val=lambda)
347 120 : CALL cite_reference(Lee1988)
348 :
349 : CALL xc_rho_set_get(rho_set, &
350 : rhoa=rhoa, rhob=rhob, norm_drhoa=norm_drhoa, &
351 : norm_drhob=norm_drhob, norm_drho=norm_drho, &
352 : rho_cutoff=epsilon_rho, &
353 120 : local_bounds=bo)
354 120 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
355 :
356 120 : dummy => rhoa
357 120 : e_0 => dummy
358 120 : e_ra => dummy
359 120 : e_rb => dummy
360 120 : e_ndra_ra => dummy
361 120 : e_ndra_rb => dummy
362 120 : e_ndrb_ra => dummy
363 120 : e_ndrb_rb => dummy
364 120 : e_ndr_ndr => dummy
365 120 : e_ndra_ndra => dummy
366 120 : e_ndrb_ndrb => dummy
367 120 : e_ndr => dummy
368 120 : e_ndra => dummy
369 120 : e_ndrb => dummy
370 120 : e_ra_ra => dummy
371 120 : e_ra_rb => dummy
372 120 : e_rb_rb => dummy
373 120 : e_ndr_ra => dummy
374 120 : e_ndr_rb => dummy
375 :
376 120 : IF (grad_deriv >= 0) THEN
377 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
378 120 : allocate_deriv=.TRUE.)
379 120 : CALL xc_derivative_get(deriv, deriv_data=e_0)
380 : END IF
381 120 : IF (grad_deriv == 1 .OR. grad_deriv == -1) THEN
382 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
383 76 : allocate_deriv=.TRUE.)
384 76 : CALL xc_derivative_get(deriv, deriv_data=e_ra)
385 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
386 76 : allocate_deriv=.TRUE.)
387 76 : CALL xc_derivative_get(deriv, deriv_data=e_rb)
388 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
389 76 : allocate_deriv=.TRUE.)
390 76 : CALL xc_derivative_get(deriv, deriv_data=e_ndr)
391 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
392 76 : allocate_deriv=.TRUE.)
393 76 : CALL xc_derivative_get(deriv, deriv_data=e_ndra)
394 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
395 76 : allocate_deriv=.TRUE.)
396 76 : CALL xc_derivative_get(deriv, deriv_data=e_ndrb)
397 : END IF
398 120 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
399 0 : CPABORT("derivatives bigger than 1 not implemented")
400 : END IF
401 :
402 : !$OMP PARALLEL DEFAULT(NONE) &
403 : !$OMP SHARED(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob) &
404 : !$OMP SHARED(e_0, e_ra, e_rb, e_ndr, e_ndra, e_ndrb) &
405 : !$OMP SHARED(grad_deriv, npoints) &
406 120 : !$OMP SHARED(epsilon_rho, lambda)
407 :
408 : CALL lyp_adiabatic_lsd_calc( &
409 : rhoa=rhoa, rhob=rhob, norm_drho=norm_drho, norm_drhoa=norm_drhoa, &
410 : norm_drhob=norm_drhob, e_0=e_0, e_ra=e_ra, e_rb=e_rb, &
411 : e_ndr=e_ndr, &
412 : e_ndra=e_ndra, e_ndrb=e_ndrb, &
413 : grad_deriv=grad_deriv, npoints=npoints, &
414 : epsilon_rho=epsilon_rho, lambda=lambda)
415 :
416 : !$OMP END PARALLEL
417 :
418 120 : CALL timestop(handle)
419 120 : END SUBROUTINE lyp_adiabatic_lsd_eval
420 :
421 : ! **************************************************************************************************
422 : !> \brief ...
423 : !> \param rhoa ...
424 : !> \param rhob ...
425 : !> \param norm_drho ...
426 : !> \param norm_drhoa ...
427 : !> \param norm_drhob ...
428 : !> \param e_0 ...
429 : !> \param e_ra ...
430 : !> \param e_rb ...
431 : !> \param e_ndr ...
432 : !> \param e_ndra ...
433 : !> \param e_ndrb ...
434 : !> \param grad_deriv ...
435 : !> \param npoints ...
436 : !> \param epsilon_rho ...
437 : !> \param lambda ...
438 : !> \par History
439 : !> 08.2008 created [mguidon]
440 : !> \author Manuel Guidon
441 : ! **************************************************************************************************
442 120 : SUBROUTINE lyp_adiabatic_lsd_calc(rhoa, rhob, norm_drho, norm_drhoa, norm_drhob, &
443 : e_0, e_ra, e_rb, &
444 : e_ndr, &
445 : e_ndra, e_ndrb, &
446 : grad_deriv, npoints, epsilon_rho, lambda)
447 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rhoa, rhob, norm_drho, norm_drhoa, &
448 : norm_drhob
449 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_ra, e_rb, e_ndr, e_ndra, e_ndrb
450 : INTEGER, INTENT(in) :: grad_deriv, npoints
451 : REAL(kind=dp), INTENT(in) :: epsilon_rho, lambda
452 :
453 : INTEGER :: ii
454 : REAL(KIND=dp) :: cf, my_ndrho, my_ndrhoa, my_ndrhob, my_rhoa, my_rhob, t1, t10, t100, t102, &
455 : t103, t106, t108, t113, t115, t118, t119, t124, t125, t128, t129, t132, t135, t138, t14, &
456 : t140, t141, t143, t145, t146, t15, t151, t153, t157, t16, t162, t165, t169, t17, t171, &
457 : t174, t179, t18, t183, t186, t187, t188, t19, t194, t196, t199, t2, t200, t202, t21, &
458 : t212, t216, t220, t222, t223, t225, t23, t231, t237, t24, t246, t25, t250, t259, t26, &
459 : t266, t27, t270, t273, t276, t28, t280, t285, t288, t294, t3, t30, t300, t307, t31, t316, &
460 : t32, t325, t348, t351, t355, t362, t387, t39, t394, t4, t41, t42
461 : REAL(KIND=dp) :: t421, t46, t47, t48, t49, t5, t51, t55, t58, t6, t62, t63, t65, t67, t7, &
462 : t73, t74, t76, t77, t78, t80, t83, t84, t85, t86, t87, t9, t90, t91, t94, t95, t96, t97
463 :
464 120 : cf = 0.3_dp*(3._dp*pi*pi)**(2._dp/3._dp)
465 :
466 120 : !$OMP DO
467 :
468 : DO ii = 1, npoints
469 4873920 : my_rhoa = MAX(rhoa(ii), 0.0_dp)
470 4873920 : my_rhob = MAX(rhob(ii), 0.0_dp)
471 4873920 : IF (my_rhoa + my_rhob > epsilon_rho) THEN
472 4862520 : my_ndrhoa = norm_drhoa(ii)
473 4862520 : my_ndrhob = norm_drhob(ii)
474 4862520 : my_ndrho = norm_drho(ii)
475 4862520 : IF (grad_deriv >= 0) THEN
476 4862520 : t1 = a*my_rhoa
477 4862520 : t2 = my_rhoa + my_rhob
478 4862520 : t3 = 0.1e1_dp/t2
479 4862520 : t4 = my_rhob*t3
480 4862520 : t5 = d*lambda
481 4862520 : t6 = t2**(0.1e1_dp/0.3e1_dp)
482 4862520 : t7 = 0.1e1_dp/t6
483 4862520 : t9 = 0.10e1_dp + t5*t7
484 4862520 : t10 = 0.1e1_dp/t9
485 4862520 : t14 = a*b
486 4862520 : t15 = c*lambda
487 4862520 : t16 = t15*t7
488 4862520 : t17 = EXP(-t16)
489 4862520 : t18 = t14*t17
490 4862520 : t19 = t2**2
491 4862520 : t21 = t6**2
492 4862520 : t23 = 0.1e1_dp/t21/t19/t2
493 4862520 : t24 = t10*t23
494 4862520 : t25 = my_rhoa*my_rhob
495 4862520 : t26 = my_rhoa**2
496 4862520 : t27 = my_rhoa**(0.1e1_dp/0.3e1_dp)
497 4862520 : t28 = t27**2
498 4862520 : t30 = my_rhob**2
499 4862520 : t31 = my_rhob**(0.1e1_dp/0.3e1_dp)
500 4862520 : t32 = t31**2
501 4862520 : t39 = t5*t7*t10
502 : t41 = 0.26111111111111111111e1_dp - 0.3888888889e0_dp*t16 - 0.3888888889e0_dp &
503 4862520 : *t39
504 4862520 : t42 = my_ndrho**2
505 : t46 = 0.25000000000000000000e1_dp - 0.5555555556e-1_dp*t16 - 0.5555555556e-1_dp &
506 4862520 : *t39
507 4862520 : t47 = my_ndrhoa**2
508 4862520 : t48 = my_ndrhob**2
509 4862520 : t49 = t47 + t48
510 4862520 : t51 = t16 + t39 - 0.110e2_dp
511 4862520 : t55 = my_rhoa*t3*t47 + t4*t48
512 : t58 = 0.12699208415745595798e2_dp*Cf*(t28*t26 + t32*t30) + t41 &
513 4862520 : *t42 - t46*t49 - 0.1111111111e0_dp*t51*t55
514 4862520 : t62 = 0.66666666666666666667e0_dp*t19
515 4862520 : t63 = t62 - t26
516 4862520 : t65 = t62 - t30
517 4862520 : t67 = t25*t58 - 0.6666666667e0_dp*t19*t42 + t63*t48 + t65*t47
518 4862520 : t73 = lambda**2
519 4862520 : t74 = t1*my_rhob
520 4862520 : t76 = 0.1e1_dp/t6/t2
521 4862520 : t77 = t9**2
522 4862520 : t78 = 0.1e1_dp/t77
523 4862520 : t80 = t76*t78*d
524 4862520 : t83 = t14*c
525 4862520 : t84 = t19**2
526 4862520 : t85 = 0.1e1_dp/t84
527 4862520 : t86 = t85*t17
528 4862520 : t87 = t10*t67
529 4862520 : t90 = t78*t85
530 4862520 : t91 = t67*d
531 4862520 : t94 = t17*t10
532 4862520 : t95 = t14*t94
533 4862520 : t96 = t23*my_rhoa
534 4862520 : t97 = c*t7
535 4862520 : t100 = d*t7*t10
536 4862520 : t102 = d**2
537 4862520 : t103 = t102*lambda
538 4862520 : t106 = t103/t21*t78
539 : t108 = -0.3888888889e0_dp*t97 - 0.3888888889e0_dp*t100 + 0.38888888888888888889e0_dp &
540 4862520 : *t106
541 : t113 = -0.5555555556e-1_dp*t97 - 0.5555555556e-1_dp*t100 + 0.55555555555555555556e-1_dp &
542 4862520 : *t106
543 4862520 : t115 = t97 + t100 - t106
544 4862520 : t118 = t108*t42 - t113*t49 - 0.1111111111e0_dp*t115*t55
545 4862520 : t119 = my_rhob*t118
546 :
547 : e_0(ii) = e_0(ii) + 0.20e1_dp*lambda*(-0.40e1_dp*t1*t4*t10 - t18*t24*t67) &
548 : + t73*(0.40e1_dp*t74*t80 + t83*t86*t87 + t18*t90*t91 - &
549 4862520 : t95*t96*t119)
550 :
551 : END IF
552 4862520 : IF (grad_deriv == 1 .OR. grad_deriv == -1) THEN
553 1398312 : t124 = a*my_rhob
554 1398312 : t125 = t3*t10
555 1398312 : t128 = 0.1e1_dp/t19
556 1398312 : t129 = my_rhob*t128
557 1398312 : t132 = 0.40e1_dp*t1*t129*t10
558 1398312 : t135 = 0.1e1_dp/t6/t19*t78
559 1398312 : t138 = 0.1333333333e1_dp*t74*t135*t5
560 1398312 : t140 = t84*t2
561 1398312 : t141 = 0.1e1_dp/t140
562 1398312 : t143 = t141*t17*t87
563 1398312 : t145 = t14*t15*t143/0.3e1_dp
564 1398312 : t146 = t17*t78
565 1398312 : t151 = t14*t146*t141*t67*t5/0.3e1_dp
566 1398312 : t153 = 0.1e1_dp/t21/t84
567 1398312 : t157 = 0.11e2_dp/0.3e1_dp*t18*t10*t153*t67
568 1398312 : t162 = t15*t76
569 1398312 : t165 = t5*t76*t10
570 1398312 : t169 = 0.1e1_dp/t21/t2
571 1398312 : t171 = t102*t73*t169*t78
572 : t174 = (0.12962962962962962963e0_dp*t162 + 0.12962962962962962963e0_dp &
573 1398312 : *t165 - 0.1296296296e0_dp*t171)*t42
574 : t179 = (0.18518518518518518519e-1_dp*t162 + 0.18518518518518518519e-1_dp &
575 1398312 : *t165 - 0.1851851852e-1_dp*t171)*t49
576 : t183 = 0.1111111111e0_dp*(-t162/0.3e1_dp - t165/0.3e1_dp + t171/0.3e1_dp) &
577 1398312 : *t55
578 1398312 : t186 = my_rhoa*t128*t47
579 1398312 : t187 = t129*t48
580 1398312 : t188 = t3*t47 - t186 - t187
581 1398312 : t194 = 0.1333333333e1_dp*t2*t42
582 1398312 : t196 = 0.13333333333333333333e1_dp*my_rhob
583 1398312 : t199 = 0.13333333333333333333e1_dp*my_rhoa
584 1398312 : t200 = t199 + t196
585 : t202 = my_rhob*t58 + t25*(0.33864555775321588795e2_dp*Cf*t28*my_rhoa &
586 : + t174 - t179 - t183 - 0.1111111111e0_dp*t51*t188) - t194 + (-0.6666666667e0_dp &
587 1398312 : *my_rhoa + t196)*t48 + t200*t47
588 1398312 : t212 = 0.5333333333e1_dp*t74*t135*d
589 1398312 : t216 = 0.1e1_dp/t77/t9
590 1398312 : t220 = 0.26666666666666666667e1_dp*t74/t21/t19*t216*t103
591 1398312 : t222 = 4*t83*t143
592 1398312 : t223 = c**2
593 1398312 : t225 = 0.1e1_dp/t6/t140
594 1398312 : t231 = t14*t223*t225*lambda*t17*t87/0.3e1_dp
595 1398312 : t237 = 0.2e1_dp/0.3e1_dp*t14*c*t225*t146*t91*lambda
596 1398312 : t246 = 0.2e1_dp/0.3e1_dp*t14*t17*t216*t225*t67*t103
597 1398312 : t250 = 4*t18*t78*t141*t91
598 1398312 : t259 = t14*t15*t141*t94*t25*t118/0.3e1_dp
599 1398312 : t266 = t14*t146*t141*t25*t118*d*lambda/0.3e1_dp
600 1398312 : t270 = 0.11e2_dp/0.3e1_dp*t95*t153*my_rhoa*t119
601 1398312 : t273 = c*t76
602 1398312 : t276 = d*t76*t10
603 1398312 : t280 = t102*t169*t78*lambda
604 1398312 : t285 = t102*d*t73*t128*t216
605 : t288 = (0.12962962962962962963e0_dp*t273 + 0.12962962962962962963e0_dp &
606 : *t276 - 0.3888888889e0_dp*t280 + 0.25925925925925925926e0_dp*t285) &
607 1398312 : *t42
608 : t294 = (0.18518518518518518519e-1_dp*t273 + 0.18518518518518518519e-1_dp &
609 : *t276 - 0.5555555556e-1_dp*t280 + 0.37037037037037037037e-1_dp*t285) &
610 1398312 : *t49
611 : t300 = 0.1111111111e0_dp*(-t273/0.3e1_dp - t276/0.3e1_dp + t280 - 0.2e1_dp &
612 1398312 : /0.3e1_dp*t285)*t55
613 : t307 = 0.40e1_dp*t124*t80 - t212 + t220 - t222 + t231 + t237 + t83 &
614 : *t86*t10*t202 + t246 - t250 + t18*t90*t202*d - t259 - &
615 : t266 + t270 - t18*t24*t119 - t95*t96*my_rhob*(t288 - t294 - &
616 1398312 : t300 - 0.1111111111e0_dp*t115*t188)
617 :
618 : e_ra(ii) = e_ra(ii) + 0.20e1_dp*lambda*(-0.40e1_dp*t124*t125 + t132 - t138 - t145 &
619 1398312 : - t151 + t157 - t18*t24*t202) + t73*t307
620 :
621 1398312 : t316 = -t186 + t3*t48 - t187
622 : t325 = my_rhoa*t58 + t25*(0.33864555775321588795e2_dp*Cf*t32*my_rhob &
623 : + t174 - t179 - t183 - 0.1111111111e0_dp*t51*t316) - t194 + t200 &
624 1398312 : *t48 + (t199 - 0.6666666667e0_dp*my_rhob)*t47
625 : t348 = 0.40e1_dp*t1*t80 - t212 + t220 - t222 + t231 + t237 + t83* &
626 : t86*t10*t325 + t246 - t250 + t18*t90*t325*d - t259 - t266 &
627 : + t270 - t18*t24*my_rhoa*t118 - t95*t96*my_rhob*(t288 - t294 &
628 1398312 : - t300 - 0.1111111111e0_dp*t115*t316)
629 :
630 : e_rb(ii) = e_rb(ii) + 0.20e1_dp*lambda*(-0.40e1_dp*t1*t125 + t132 - t138 - t145 - &
631 1398312 : t151 + t157 - t18*t24*t325) + t73*t348
632 :
633 1398312 : t351 = lambda*a*b
634 1398312 : t355 = t3*my_ndrhoa
635 : t362 = t25*(-REAL(2*t46*my_ndrhoa, dp) - 0.2222222222e0_dp*t51*my_rhoa &
636 1398312 : *t355) + REAL(2*t65*my_ndrhoa, dp)
637 :
638 : e_ndra(ii) = e_ndra(ii) - 0.20e1_dp*t351*t94*t23*t362 + t73*(t83*t86*t10* &
639 : t362 + t18*t90*t362*d - t95*t96*my_rhob*(-REAL(2*t113* &
640 1398312 : my_ndrhoa, dp) - 0.2222222222e0_dp*t115*my_rhoa*t355))
641 :
642 1398312 : t387 = t3*my_ndrhob
643 : t394 = t25*(-REAL(2*t46*my_ndrhob, dp) - 0.2222222222e0_dp*t51*my_rhob &
644 1398312 : *t387) + REAL(2*t63*my_ndrhob, dp)
645 :
646 : e_ndrb(ii) = e_ndrb(ii) - 0.20e1_dp*t351*t94*t23*t394 + t73*(t83*t86*t10* &
647 : t394 + t18*t90*t394*d - t95*t96*my_rhob*(-REAL(2*t113* &
648 1398312 : my_ndrhob, dp) - 0.2222222222e0_dp*t115*my_rhob*t387))
649 :
650 1398312 : t421 = REAL(2*t25*t41*my_ndrho, dp) - 0.1333333333e1_dp*REAL(t19, dp)*REAL(my_ndrho, dp)
651 :
652 : e_ndr(ii) = e_ndr(ii) - 0.20e1_dp*t351*t94*t23*t421 + t73*(t83*t86*t10*t421 &
653 1398312 : + t18*t90*t421*d - REAL(2*t95*t96*my_rhob*t108*my_ndrho, dp))
654 :
655 : END IF
656 : END IF
657 : END DO
658 :
659 : !$OMP END DO
660 :
661 120 : END SUBROUTINE lyp_adiabatic_lsd_calc
662 :
663 : END MODULE xc_lyp_adiabatic
|
