(git:34ef472)
xc_tpss.F
Go to the documentation of this file.
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 tpss functional.
10 !> \note
11 !> The derivation of the formulaes is lengthly, and not fully trivial,
12 !> so I have put it in doc/tpss.mw
13 !> \par History
14 !> 05.2004 created
15 !> \author fawzi
16 ! **************************************************************************************************
17 MODULE xc_tpss
18  USE bibliography, ONLY: tao2003,&
19  cite_reference
20  USE cp_log_handling, ONLY: cp_get_default_logger,&
21  cp_logger_type
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_rho,&
28  deriv_tau
29  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
30  xc_dset_get_derivative
31  USE xc_derivative_types, ONLY: xc_derivative_get,&
32  xc_derivative_type
33  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
34  USE xc_rho_set_types, ONLY: xc_rho_set_get,&
35  xc_rho_set_type
36 #include "../base/base_uses.f90"
37 
38  IMPLICIT NONE
39  PRIVATE
40 
41  LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .true.
42  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_tpss'
43 
44  PUBLIC :: tpss_lda_info, tpss_lda_eval
45 
46 !***
47 CONTAINS
48 
49 ! **************************************************************************************************
50 !> \brief return various information on the functional
51 !> \param tpss_params ...
52 !> \param reference string with the reference of the actual functional
53 !> \param shortform string with the shortform of the functional name
54 !> \param needs the components needed by this functional are set to
55 !> true (does not set the unneeded components to false)
56 !> \param max_deriv the highest derivative available
57 !> \author fawzi
58 ! **************************************************************************************************
59  SUBROUTINE tpss_lda_info(tpss_params, reference, shortform, needs, max_deriv)
60  TYPE(section_vals_type), POINTER :: tpss_params
61  CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
62  TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
63  INTEGER, INTENT(out), OPTIONAL :: max_deriv
64 
65  REAL(kind=dp) :: sc, sx
66 
67  CALL section_vals_val_get(tpss_params, "SCALE_C", r_val=sc)
68  CALL section_vals_val_get(tpss_params, "SCALE_X", r_val=sx)
69 
70  IF (PRESENT(reference)) THEN
71  IF (sx == 1._dp .AND. sc == 1._dp) THEN
72  reference = "J. Tao, J.P.Perdew, V.N.Staroverov, E.Scuseria PRL, 91, 146401 (2003) {LDA version}"
73  ELSE
74  WRITE (reference, "(a,'sx=',f5.3,'sc=',f5.3,' {LDA version}')") &
75  "J. Tao, J.P.Perdew, V.N.Staroverov, E.Scuseria PRL, 91, 146401 (2003)", &
76  sx, sc
77  END IF
78  END IF
79  IF (PRESENT(shortform)) THEN
80  IF (sx == 1._dp .AND. sc == 1._dp) THEN
81  shortform = "TPSS meta-GGA functional (LDA)"
82  ELSE
83  WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,' (LDA)')") &
84  "TPSS meta-GGA functional", &
85  sx, sc
86  END IF
87  END IF
88  IF (PRESENT(needs)) THEN
89  needs%rho = .true.
90  needs%tau = .true.
91  needs%norm_drho = .true.
92  END IF
93  IF (PRESENT(max_deriv)) max_deriv = 1
94 
95  END SUBROUTINE tpss_lda_info
96 
97 ! **************************************************************************************************
98 !> \brief evaluates the tpss functional in the spin unpolarized (lda) case
99 !> \param rho_set the density where you want to evaluate the functional
100 !> \param deriv_set place where to store the functional derivatives (they are
101 !> added to the derivatives)
102 !> \param grad_deriv degree of the derivative that should be evaluated,
103 !> if positive all the derivatives up to the given degree are evaluated,
104 !> if negative only the given degree is calculated
105 !> \param tpss_params ...
106 !> \author fawzi
107 ! **************************************************************************************************
108  SUBROUTINE tpss_lda_eval(rho_set, deriv_set, grad_deriv, tpss_params)
109  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
110  TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
111  INTEGER, INTENT(in) :: grad_deriv
112  TYPE(section_vals_type), POINTER :: tpss_params
113 
114  CHARACTER(len=*), PARAMETER :: routinen = 'tpss_lda_eval'
115 
116  INTEGER :: handle, non_coer, npoints
117  INTEGER, DIMENSION(2, 3) :: bo
118  REAL(kind=dp) :: epsilon_rho, epsilon_tau, scale_ec, &
119  scale_ex
120  REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
121  POINTER :: dummy, e_0, e_ndrho, e_rho, e_tau, &
122  norm_drho, rho, tau
123  TYPE(cp_logger_type), POINTER :: logger
124  TYPE(xc_derivative_type), POINTER :: deriv
125 
126  CALL timeset(routinen, handle)
127 
128  CALL cite_reference(tao2003)
129 
130  CALL xc_rho_set_get(rho_set, rho=rho, &
131  norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
132  tau=tau, tau_cutoff=epsilon_tau)
133  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
134 
135  dummy => rho
136 
137  e_0 => dummy
138  e_rho => dummy
139  e_ndrho => dummy
140  e_tau => dummy
141 
142  IF (grad_deriv >= 0) THEN
143  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
144  allocate_deriv=.true.)
145  CALL xc_derivative_get(deriv, deriv_data=e_0)
146  END IF
147  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
148  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
149  allocate_deriv=.true.)
150  CALL xc_derivative_get(deriv, deriv_data=e_rho)
151  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
152  allocate_deriv=.true.)
153  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
154  deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
155  allocate_deriv=.true.)
156  CALL xc_derivative_get(deriv, deriv_data=e_tau)
157  END IF
158  IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
159  cpabort("derivatives bigger than 1 not implemented")
160  END IF
161 
162  non_coer = 0
163  CALL section_vals_val_get(tpss_params, "SCALE_C", r_val=scale_ec)
164  CALL section_vals_val_get(tpss_params, "SCALE_X", r_val=scale_ex)
165 
166 !$OMP PARALLEL DEFAULT(NONE) &
167 !$OMP SHARED(rho, tau, norm_drho, e_0, e_rho, e_ndrho, e_tau) &
168 !$OMP SHARED(epsilon_rho, epsilon_tau, npoints, grad_deriv) &
169 !$OMP SHARED(scale_ec, scale_ex) &
170 !$OMP REDUCTION(+: non_coer)
171 
172  CALL tpss_lda_calc(rho=rho, norm_drho=norm_drho, &
173  tau=tau, e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_tau=e_tau, &
174  grad_deriv=grad_deriv, npoints=npoints, epsilon_rho=epsilon_rho, &
175  epsilon_tau=epsilon_tau, scale_ec=scale_ec, scale_ex=scale_ex, non_coer=non_coer)
176 
177 !$OMP END PARALLEL
178 
179  logger => cp_get_default_logger()
180  ! we could check if tau/grad were consistent, but don't do anything here
181  IF (non_coer > 0) THEN
182  non_coer = 0
183  END IF
184 
185  CALL timestop(handle)
186  END SUBROUTINE tpss_lda_eval
187 
188 ! **************************************************************************************************
189 !> \brief low level calculation routine for the unpolarized (lda) tpss
190 !> \param rho ...
191 !> \param norm_drho ...
192 !> \param tau ...
193 !> \param e_0 ...
194 !> \param e_rho ...
195 !> \param e_ndrho ...
196 !> \param e_tau ...
197 !> \param npoints ...
198 !> \param grad_deriv ...
199 !> \param epsilon_rho ...
200 !> \param epsilon_tau ...
201 !> \param scale_ec ...
202 !> \param scale_ex ...
203 !> \param non_coer ...
204 !> \author fawzi
205 !> \note
206 !> maple is nice, but if you want the uman readable version of the code
207 !> look in doc/tpss.mw
208 ! **************************************************************************************************
209  SUBROUTINE tpss_lda_calc(rho, norm_drho, tau, e_0, e_rho, e_ndrho, e_tau, &
210  npoints, grad_deriv, epsilon_rho, epsilon_tau, &
211  scale_ec, scale_ex, non_coer)
212  REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho, norm_drho, tau
213  REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho, e_ndrho, e_tau
214  INTEGER, INTENT(in) :: npoints, grad_deriv
215  REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_tau, scale_ec, &
216  scale_ex
217  INTEGER, INTENT(inout) :: non_coer
218 
219  INTEGER :: abs_grad_deriv, ii
220  LOGICAL :: t571, t639
221  REAL(kind=dp) :: a, a_1, a_2, a_s1, a_s1rho, a_s2, a_s2rho, alpha, alpha_1_1, alpha_1_2, &
222  alphanorm_drho, alpharho, alphatau, arho, b, beta, beta_1_1, beta_1_2, beta_2_1, &
223  beta_2_2, beta_3_1, beta_3_2, beta_4_1, beta_4_2, c, d, e_c_u_0, e_c_u_0rho, e_c_u_1_s1, &
224  e_c_u_1_s1rho, e_c_u_1_s2, e_c_u_1_s2rho, e_var, epsilon_cgga, epsilon_cgga_0_1, &
225  epsilon_cgga_1_0, epsilon_cggarho, epsilon_crevpkzb, epsilon_crevpkzbnorm_drho, &
226  epsilon_crevpkzbrho, epsilon_crevpkzbtau, ex_unif, fx, gamma_var, hnorm_drho, k_f_s1, &
227  k_f_s1rho, k_s, k_s_s1, k_s_s2, kappa, m, ma, manorm_drho, marho, mb, mbnorm_drho, mbrho
228  REAL(kind=dp) :: mu, my_ndrho, my_rho, my_tau, p, p_1, p_2, p_3, phi_s1, phi_s2, pnorm_drho, &
229  prho, rs, rs_s1, rs_s1rho, rs_s2, rs_s2rho, rsrho, t, t1, t100, t101, t111, t12, t13, &
230  t138, t14, t140, t143, t145, t146, t147, t151, t152, t16, t161, t168, t177, t186, t187, &
231  t189, t19, t190, t191, t193, t194, t196, t197, t198, t199, t2, t20, t201, t202, t204, &
232  t205, t208, t209, t21, t211, t212, t213, t215, t216, t218, t219, t22, t220, t221, t223, &
233  t224, t226, t227, t230, t231, t233, t234, t235, t238, t239, t241, t242, t243, t245, t246, &
234  t248, t249, t252, t253, t254, t256, t26, t260, t263, t264, t265
235  REAL(kind=dp) :: t267, t268, t269, t27, t271, t272, t274, t275, t276, t277, t278, t279, t28, &
236  t280, t281, t284, t286, t288, t29, t290, t291, t293, t294, t295, t299, t3, t301, t302, &
237  t303, t305, t307, t310, t313, t316, t319, t322, t325, t327, t328, t329, t331, t337, t340, &
238  t343, t344, t35, t351, t36, t370, t371, t376, t383, t385, t386, t39, t390, t391, t395, &
239  t396, t398, t4, t403, t404, t406, t41, t410, t411, t419, t42, t430, t437, t445, t450, &
240  t452, t464, t472, t475, t485, t489, t49, t490, t5, t505, t513, t517, t536, t541, t542, &
241  t546, t547, t549, t55, t554, t555, t557, t561, t562, t569, t574
242  REAL(kind=dp) :: t58, t585, t6, t60, t604, t609, t610, t614, t615, t617, t622, t623, t625, &
243  t629, t630, t637, t642, t645, t659, t67, t7, t71, t73, t77, t78, t79, t799, t80, t84, &
244  t85, t89, t9, t94, t95, t96, t_s1, t_s1norm_drho, t_s1rho, t_s2, t_s2norm_drho, t_s2rho, &
245  tau_w, tau_wnorm_drho, tau_wrho, tildeq_b, tildeq_bnorm_drho, tildeq_brho, tildeq_btau, &
246  tnorm_drho, trho, z, znorm_drho, zrho, ztau
247 
248  IF (.false.) THEN
249  ! useful for testing, we just hack in a well defined functional of tau, ndrho and rho
250  ! and see that things converge properly with OT.
251 !$OMP DO
252  DO ii = 1, npoints
253  my_tau = tau(ii)
254  my_rho = rho(ii)
255  my_ndrho = norm_drho(ii)
256  IF (grad_deriv >= 0) THEN
257  e_0(ii) = e_0(ii) + my_tau*my_ndrho*my_rho
258  END IF
259  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
260  e_rho(ii) = e_rho(ii) + my_tau*my_ndrho
261  e_ndrho(ii) = e_ndrho(ii) + my_tau*my_rho
262  e_tau(ii) = e_tau(ii) + my_rho*my_ndrho
263  END IF
264  END DO
265 !$OMP END DO
266  RETURN
267  END IF
268 
269  abs_grad_deriv = abs(grad_deriv)
270 
271  kappa = 0.804e0_dp
272  beta = 0.66725e-1_dp
273  mu = 0.21951e0_dp
274  gamma_var = (0.1e1_dp - log(0.2e1_dp))/pi**2
275  b = 0.4e0_dp
276  c = 0.159096e1_dp
277  e_var = 0.1537e1_dp
278  d = 0.28e1_dp
279  p_1 = 0.10e1_dp
280  a_1 = 0.31091e-1_dp
281  alpha_1_1 = 0.21370e0_dp
282  beta_1_1 = 0.75957e1_dp
283  beta_2_1 = 0.35876e1_dp
284  beta_3_1 = 0.16382e1_dp
285  beta_4_1 = 0.49294e0_dp
286  p_2 = 0.10e1_dp
287  a_2 = 0.15545e-1_dp
288  alpha_1_2 = 0.20548e0_dp
289  beta_1_2 = 0.141189e2_dp
290  beta_2_2 = 0.61977e1_dp
291  beta_3_2 = 0.33662e1_dp
292  beta_4_2 = 0.62517e0_dp
293  p_3 = 0.10e1_dp
294 
295  t1 = 3._dp**(0.1e1_dp/0.3e1_dp)
296  t2 = 4._dp**(0.1e1_dp/0.3e1_dp)
297  t3 = t2**2
298  t4 = t1*t3
299  t5 = 2._dp**(0.1e1_dp/0.3e1_dp)
300  t6 = 0.1e1_dp/pi
301  t12 = t5**2
302 
303 !$OMP DO
304 
305  DO ii = 1, npoints
306  my_tau = tau(ii)
307  my_rho = rho(ii)
308  IF (my_rho > epsilon_rho .AND. my_tau > epsilon_tau) THEN
309  my_ndrho = norm_drho(ii)
310 
311  t7 = 0.1e1_dp/my_rho
312  t254 = my_ndrho**2
313  tau_w = t254*t7/0.8e1_dp
314 
315  IF (my_tau < tau_w) THEN
316  ! enforce z=norm_rho**2/(8._dp*rho*tau) <1
317  m = 0.5_dp*t254 + 4.0_dp*my_rho*my_tau
318  my_tau = m/8._dp/my_rho
319  my_ndrho = sqrt(m)
320  t254 = m
321  non_coer = non_coer + 1
322  END IF
323 
324  t9 = (t6*t7)**(0.1e1_dp/0.3e1_dp)
325  rs_s1 = t4*t5*t9/0.4e1_dp
326  phi_s1 = t12/0.2e1_dp
327  t13 = t1*t12
328  t14 = pi**2
329  t16 = (t14*my_rho)**(0.1e1_dp/0.3e1_dp)
330  k_f_s1 = t13*t16/0.2e1_dp
331  t19 = sqrt(k_f_s1*t6)
332  k_s_s1 = 0.2e1_dp*t19
333  t20 = 0.1e1_dp/phi_s1
334  t21 = my_ndrho*t20
335  t22 = 0.1e1_dp/k_s_s1
336  t_s1 = t21*t22*t7/0.2e1_dp
337  rs_s2 = rs_s1
338  phi_s2 = phi_s1
339  t26 = sqrt(k_f_s1*t6)
340  k_s_s2 = 0.2e1_dp*t26
341  t27 = 0.1e1_dp/phi_s2
342  t28 = my_ndrho*t27
343  t29 = 0.1e1_dp/k_s_s2
344  t_s2 = t28*t29*t7/0.2e1_dp
345  t35 = 0.1e1_dp/a_1
346  t36 = sqrt(rs_s2)
347  t39 = t36*rs_s2
348  t41 = p_1 + 0.1e1_dp
349  t42 = rs_s2**t41
350  t49 = log(0.1e1_dp + t35/(beta_1_1*t36 + beta_2_1*rs_s2 + &
351  beta_3_1*t39 + beta_4_1*t42)/0.2e1_dp)
352  t55 = sqrt(rs_s1)
353  t58 = t55*rs_s1
354  t60 = rs_s1**t41
355  t67 = log(0.1e1_dp + t35/(beta_1_1*t55 + beta_2_1*rs_s1 + &
356  beta_3_1*t58 + beta_4_1*t60)/0.2e1_dp)
357  t71 = 0.1e1_dp + alpha_1_2*rs_s2
358  t73 = 0.1e1_dp/a_2
359  t77 = p_2 + 0.1e1_dp
360  t78 = rs_s2**t77
361  t79 = beta_4_2*t78
362  t80 = beta_1_2*t36 + beta_2_2*rs_s2 + beta_3_2*t39 + t79
363  t84 = 0.1e1_dp + t73/t80/0.2e1_dp
364  t85 = log(t84)
365  e_c_u_1_s2 = -0.2e1_dp*a_2*t71*t85
366  t89 = 0.1e1_dp + alpha_1_2*rs_s1
367  t94 = rs_s1**t77
368  t95 = beta_4_2*t94
369  t96 = beta_1_2*t55 + beta_2_2*rs_s1 + beta_3_2*t58 + t95
370  t100 = 0.1e1_dp + t73/t96/0.2e1_dp
371  t101 = log(t100)
372  e_c_u_1_s1 = -0.2e1_dp*a_2*t89*t101
373  t111 = p_3 + 1._dp
374  rs = t4*t9/0.4e1_dp
375  t138 = 0.1e1_dp + alpha_1_1*rs
376  t140 = sqrt(rs)
377  t143 = t140*rs
378  t145 = rs**t41
379  t146 = beta_4_1*t145
380  t147 = beta_1_1*t140 + beta_2_1*rs + beta_3_1*t143 + t146
381  t151 = 0.1e1_dp + t35/t147/0.2e1_dp
382  t152 = log(t151)
383  e_c_u_0 = -0.2e1_dp*a_1*t138*t152
384  t161 = rs**t77
385  t168 = log(0.1e1_dp + t73/(beta_1_2*t140 + beta_2_2*rs + &
386  beta_3_2*t143 + beta_4_2*t161)/0.2e1_dp)
387  t177 = rs**t111
388  t186 = 0.1e1_dp/gamma_var
389  t187 = beta*t186
390  t189 = phi_s1**2
391  t190 = t189*phi_s1
392  t191 = 0.1e1_dp/t190
393  t193 = exp(-e_c_u_1_s1*t186*t191)
394  t194 = t193 - 0.1e1_dp
395  a_s1 = t187/t194
396  t196 = gamma_var*t190
397  t197 = t_s1**2
398  t198 = a_s1*t197
399  t199 = 0.1e1_dp + t198
400  t201 = a_s1**2
401  t202 = t197**2
402  t204 = 0.1e1_dp + t198 + t201*t202
403  t205 = 0.1e1_dp/t204
404  t208 = 0.1e1_dp + t187*t197*t199*t205
405  t209 = log(t208)
406  epsilon_cgga_1_0 = e_c_u_1_s1 + t196*t209
407  t211 = phi_s2**2
408  t212 = t211*phi_s2
409  t213 = 0.1e1_dp/t212
410  t215 = exp(-e_c_u_1_s2*t186*t213)
411  t216 = t215 - 0.1e1_dp
412  a_s2 = t187/t216
413  t218 = gamma_var*t212
414  t219 = t_s2**2
415  t220 = a_s2*t219
416  t221 = t220 + 0.1e1_dp
417  t223 = a_s2**2
418  t224 = t219**2
419  t226 = 0.1e1_dp + t220 + t223*t224
420  t227 = 0.1e1_dp/t226
421  t230 = 0.1e1_dp + t187*t219*t221*t227
422  t231 = log(t230)
423  epsilon_cgga_0_1 = e_c_u_1_s2 + t218*t231
424  t233 = sqrt(t1*t16*t6)
425  k_s = 0.2e1_dp*t233
426  t234 = 0.1e1_dp/k_s
427  t235 = my_ndrho*t234
428  t = t235*t7/0.2e1_dp
429  t238 = exp(-e_c_u_0*t186)
430  t239 = -0.1e1_dp + t238
431  a = t187/t239
432  t241 = t**2
433  t242 = a*t241
434  t243 = 0.1e1_dp + t242
435  t245 = a**2
436  t246 = t241**2
437  t248 = 0.1e1_dp + t242 + t245*t246
438  t249 = 0.1e1_dp/t248
439  t252 = 0.1e1_dp + t187*t241*t243*t249
440  t253 = log(t252)
441  epsilon_cgga = e_c_u_0 + gamma_var*t253
442  ma = max(epsilon_cgga_1_0, epsilon_cgga)
443  mb = max(epsilon_cgga_0_1, epsilon_cgga)
444  t256 = tau_w**2
445  t260 = ma/0.2e1_dp + mb/0.2e1_dp
446  t263 = 0.53e0_dp*epsilon_cgga*t256 - 0.153e1_dp*t256*t260
447  t264 = my_tau**2
448  t265 = 0.1e1_dp/t264
449  epsilon_crevpkzb = epsilon_cgga + t263*t265
450  t267 = my_rho*epsilon_crevpkzb
451  t268 = d*epsilon_crevpkzb
452  t269 = t256*tau_w
453  t271 = 0.1e1_dp/t264/my_tau
454  t272 = t269*t271
455  t274 = 0.1e1_dp + t268*t272
456  t275 = t254*t1
457  t276 = t14**(0.1e1_dp/0.3e1_dp)
458  t277 = t276**2
459  t278 = 0.1e1_dp/t277
460  t279 = my_rho**2
461  t280 = my_rho**(0.1e1_dp/0.3e1_dp)
462  t281 = t280**2
463  t284 = t278/t281/t279
464  p = t275*t284/0.12e2_dp
465  t286 = 0.1e1_dp/my_tau
466  z = tau_w*t286
467  t288 = 0.1e1_dp/z - 0.1e1_dp
468  alpha = 0.5e1_dp/0.3e1_dp*p*t288
469  t290 = alpha - 0.1e1_dp
470  t291 = b*alpha
471  t293 = 0.1e1_dp + t291*t290
472  t294 = sqrt(t293)
473  t295 = 0.1e1_dp/t294
474  tildeq_b = 0.9e1_dp/0.20e2_dp*t290*t295 + 0.2e1_dp/0.3e1_dp*p
475  t299 = z**2
476  t301 = 0.1e1_dp + t299
477  t302 = t301**2
478  t303 = 0.1e1_dp/t302
479  t305 = 0.10e2_dp/0.81e2_dp + c*t299*t303
480  t307 = tildeq_b**2
481  t310 = p**2
482  t313 = sqrt(0.18e2_dp*t299 + 0.50e2_dp*t310)
483  t316 = 0.1e1_dp/kappa
484  t319 = sqrt(e_var)
485  t322 = e_var*mu
486  t325 = t305*p + 0.146e3_dp/0.2025e4_dp*t307 - 0.73e2_dp/ &
487  0.4050e4_dp*tildeq_b*t313 + 0.100e3_dp/0.6561e4_dp*t316* &
488  t310 + 0.4e1_dp/0.45e2_dp*t319*t299 + t322*t310*p
489  t327 = 0.1e1_dp + t319*p
490  t328 = t327**2
491  t329 = 0.1e1_dp/t328
492  t331 = 0.1e1_dp + t325*t329*t316
493  fx = 0.1e1_dp + kappa - kappa/t331
494  ex_unif = -0.3e1_dp/0.4e1_dp*t1*t16*t6
495  t337 = my_rho*ex_unif
496 
497  IF (grad_deriv >= 0) THEN
498  e_0(ii) = e_0(ii) + &
499  scale_ec*t267*t274 + scale_ex*t337*fx
500  END IF
501 
502  IF (abs_grad_deriv > 0) THEN
503  t340 = t9**2
504  t343 = 0.1e1_dp/t279
505  t344 = 0.1e1_dp/t340*t6*t343
506  rsrho = -t4*t344/0.12e2_dp
507  t351 = t147**2
508  e_c_u_0rho = -0.2e1_dp*a_1*alpha_1_1*rsrho*t152 + t138/ &
509  t351*(beta_1_1/t140*rsrho/0.2e1_dp + beta_2_1*rsrho + &
510  0.3e1_dp/0.2e1_dp*beta_3_1*t140*rsrho + t146*t41*rsrho/ &
511  rs)/t151
512  t370 = t16**2
513  t371 = 0.1e1_dp/t370
514  t376 = k_s**2
515  trho = -my_ndrho/t376*t7/t233*t1*t371*t14*t6 &
516  /0.6e1_dp - t235*t343/0.2e1_dp
517  t383 = gamma_var**2
518  t385 = beta/t383
519  t386 = t239**2
520  arho = t385/t386*e_c_u_0rho*t238
521  t390 = t187*t
522  t391 = t243*t249
523  t395 = arho*t241
524  t396 = a*t
525  t398 = 0.2e1_dp*t396*trho
526  t403 = t187*t241
527  t404 = t248**2
528  t406 = t243/t404
529  t410 = t241*t
530  t411 = t245*t410
531  t419 = 0.1e1_dp/t252
532  epsilon_cggarho = e_c_u_0rho + gamma_var*(0.2e1_dp*t390*t391 &
533  *trho + t187*t241*(t395 + t398)*t249 - t403*t406*(t395 + &
534  t398 + 0.2e1_dp*a*t246*arho + 0.4e1_dp*t411*trho))*t419
535  tau_wrho = -t254*t343/0.8e1_dp
536  prho = -0.2e1_dp/0.9e1_dp*t275*t278/t281/t279/my_rho
537  zrho = tau_wrho*t286
538  t430 = p/t299
539  alpharho = 0.5e1_dp/0.3e1_dp*prho*t288 - 0.5e1_dp/0.3e1_dp &
540  *t430*zrho
541  t437 = t290/t294/t293
542  tildeq_brho = 0.9e1_dp/0.20e2_dp*alpharho*t295 - 0.9e1_dp/ &
543  0.40e2_dp*t437*(b*alpharho*t290 + t291*alpharho) + &
544  0.2e1_dp/0.3e1_dp*prho
545  t445 = c*z
546  t450 = c*t299*z
547  t452 = 0.1e1_dp/t302/t301
548  t464 = tildeq_b/t313
549  t472 = t316*p
550  t475 = t319*z
551  t485 = t325/t328/t327
552  t489 = t331**2
553  t490 = 0.1e1_dp/t489
554  rs_s1rho = -t4*t5*t344/0.12e2_dp
555  k_f_s1rho = t13*t371*t14/0.6e1_dp
556  t505 = k_s_s1**2
557  t_s1rho = -t21/t505*t7/t19*k_f_s1rho*t6/0.2e1_dp - t21 &
558  *t22*t343/0.2e1_dp
559  t513 = a_2*alpha_1_2
560  t517 = t96**2
561  e_c_u_1_s1rho = -0.2e1_dp*t513*rs_s1rho*t101 + t89/t517* &
562  (beta_1_2/t55*rs_s1rho/0.2e1_dp + beta_2_2*rs_s1rho + &
563  0.3e1_dp/0.2e1_dp*beta_3_2*t55*rs_s1rho + t95*t77* &
564  rs_s1rho/rs_s1)/t100
565  t536 = t194**2
566  a_s1rho = t385/t536*e_c_u_1_s1rho*t191*t193
567  t541 = t187*t_s1
568  t542 = t199*t205
569  t546 = a_s1rho*t197
570  t547 = a_s1*t_s1
571  t549 = 0.2e1_dp*t547*t_s1rho
572  t554 = t187*t197
573  t555 = t204**2
574  t557 = t199/t555
575  t561 = t197*t_s1
576  t562 = t201*t561
577  t569 = 0.1e1_dp/t208
578  t571 = epsilon_cgga .LT. epsilon_cgga_1_0
579  IF (t571) THEN
580  marho = e_c_u_1_s1rho + t196*(0.2e1_dp*t541*t542 &
581  *t_s1rho + t187*t197*(t546 + t549)*t205 - t554*t557*(t546 &
582  + t549 + 0.2e1_dp*a_s1*t202*a_s1rho + 0.4e1_dp*t562* &
583  t_s1rho))*t569
584  ELSE
585  marho = epsilon_cggarho
586  END IF
587  rs_s2rho = rs_s1rho
588  t574 = k_s_s2**2
589  t_s2rho = -t28/t574*t7/t26*k_f_s1rho*t6/0.2e1_dp - t28 &
590  *t29*t343/0.2e1_dp
591  t585 = t80**2
592  e_c_u_1_s2rho = -0.2e1_dp*t513*rs_s2rho*t85 + t71/t585*( &
593  beta_1_2/t36*rs_s2rho/0.2e1_dp + beta_2_2*rs_s2rho + &
594  0.3e1_dp/0.2e1_dp*beta_3_2*t36*rs_s2rho + t79*t77* &
595  rs_s2rho/rs_s2)/t84
596  t604 = t216**2
597  a_s2rho = t385/t604*e_c_u_1_s2rho*t213*t215
598  t609 = t187*t_s2
599  t610 = t221*t227
600  t614 = a_s2rho*t219
601  t615 = a_s2*t_s2
602  t617 = 0.2e1_dp*t615*t_s2rho
603  t622 = t187*t219
604  t623 = t226**2
605  t625 = t221/t623
606  t629 = t219*t_s2
607  t630 = t223*t629
608  t637 = 0.1e1_dp/t230
609  t639 = epsilon_cgga .LT. epsilon_cgga_0_1
610  IF (t639) THEN
611  mbrho = e_c_u_1_s2rho + t218*(0.2e1_dp*t609*t610 &
612  *t_s2rho + t187*t219*(t614 + t617)*t227 - t622*t625*(t614 &
613  + t617 + 0.2e1_dp*a_s2*t224*a_s2rho + 0.4e1_dp*t630* &
614  t_s2rho))*t637
615  ELSE
616  mbrho = epsilon_cggarho
617  END IF
618  t642 = epsilon_cgga*tau_w
619  t645 = tau_w*t260
620  epsilon_crevpkzbrho = epsilon_cggarho + (0.53e0_dp* &
621  epsilon_cggarho*t256 + 0.106e1_dp*t642*tau_wrho - 0.306e1_dp* &
622  t645*tau_wrho - 0.153e1_dp*t256*(marho/0.2e1_dp + mbrho/ &
623  0.2e1_dp))*t265
624  t659 = t256*t271
625 
626  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
627  e_rho(ii) = e_rho(ii) + &
628  scale_ec*(epsilon_crevpkzb*t274 + my_rho* &
629  epsilon_crevpkzbrho*t274 + t267*(d*epsilon_crevpkzbrho*t272 &
630  + 0.3e1_dp*t268*t659*tau_wrho)) + scale_ex*(ex_unif*fx - &
631  my_rho*pi*t1*t371*fx/0.4e1_dp + t337* &
632  t490*(((0.2e1_dp*t445*t303*zrho - 0.4e1_dp*t450*t452* &
633  zrho)*p + t305*prho + 0.292e3_dp/0.2025e4_dp*tildeq_b* &
634  tildeq_brho - 0.73e2_dp/0.4050e4_dp*tildeq_brho*t313 - &
635  0.73e2_dp/0.8100e4_dp*t464*(0.36e2_dp*z*zrho + 0.100e3_dp &
636  *p*prho) + 0.200e3_dp/0.6561e4_dp*t472*prho + 0.8e1_dp/ &
637  0.45e2_dp*t475*zrho + 0.3e1_dp*t322*t310*prho)*t329 - &
638  0.2e1_dp*t485*t319*prho))
639  END IF
640 
641  tnorm_drho = t234*t7/0.2e1_dp
642  hnorm_drho = gamma_var*(0.2e1_dp*t390*t391*tnorm_drho + &
643  0.2e1_dp*t187*t410*a*tnorm_drho*t249 - t403*t406*( &
644  0.2e1_dp*t396*tnorm_drho + 0.4e1_dp*t411*tnorm_drho))*t419
645  tau_wnorm_drho = my_ndrho*t7/0.4e1_dp
646  pnorm_drho = my_ndrho*t1*t284/0.6e1_dp
647  znorm_drho = tau_wnorm_drho*t286
648  alphanorm_drho = 0.5e1_dp/0.3e1_dp*pnorm_drho*t288 - &
649  0.5e1_dp/0.3e1_dp*t430*znorm_drho
650  tildeq_bnorm_drho = 0.9e1_dp/0.20e2_dp*alphanorm_drho*t295 - &
651  0.9e1_dp/0.40e2_dp*t437*(b*alphanorm_drho*t290 + t291* &
652  alphanorm_drho) + 0.2e1_dp/0.3e1_dp*pnorm_drho
653  t_s1norm_drho = t20*t22*t7/0.2e1_dp
654  IF (t571) THEN
655  manorm_drho = t196*(0.2e1_dp*t541*t542* &
656  t_s1norm_drho + 0.2e1_dp*t187*t561*a_s1*t_s1norm_drho*t205 &
657  - t554*t557*(0.2e1_dp*t547*t_s1norm_drho + 0.4e1_dp*t562 &
658  *t_s1norm_drho))*t569
659  ELSE
660  manorm_drho = hnorm_drho
661  END IF
662  t_s2norm_drho = t27*t29*t7/0.2e1_dp
663  IF (t639) THEN
664  mbnorm_drho = t218*(0.2e1_dp*t609*t610* &
665  t_s2norm_drho + 0.2e1_dp*t187*t629*a_s2*t_s2norm_drho*t227 &
666  - t622*t625*(0.2e1_dp*t615*t_s2norm_drho + 0.4e1_dp*t630 &
667  *t_s2norm_drho))*t637
668  ELSE
669  mbnorm_drho = hnorm_drho
670  END IF
671  epsilon_crevpkzbnorm_drho = hnorm_drho + (0.53e0_dp*hnorm_drho* &
672  t256 + 0.106e1_dp*t642*tau_wnorm_drho - 0.306e1_dp*t645* &
673  tau_wnorm_drho - 0.153e1_dp*t256*(manorm_drho/0.2e1_dp + &
674  mbnorm_drho/0.2e1_dp))*t265
675 
676  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
677  e_ndrho(ii) = e_ndrho(ii) + &
678  scale_ec*(my_rho*epsilon_crevpkzbnorm_drho* &
679  t274 + t267*(d*epsilon_crevpkzbnorm_drho*t272 + 0.3e1_dp* &
680  t268*t659*tau_wnorm_drho)) + scale_ex*t337*t490*((( &
681  0.2e1_dp*t445*t303*znorm_drho - 0.4e1_dp*t450*t452* &
682  znorm_drho)*p + t305*pnorm_drho + 0.292e3_dp/0.2025e4_dp* &
683  tildeq_b*tildeq_bnorm_drho - 0.73e2_dp/0.4050e4_dp* &
684  tildeq_bnorm_drho*t313 - 0.73e2_dp/0.8100e4_dp*t464*( &
685  0.36e2_dp*z*znorm_drho + 0.100e3_dp*p*pnorm_drho) + &
686  0.200e3_dp/0.6561e4_dp*t472*pnorm_drho + 0.8e1_dp/0.45e2_dp &
687  *t475*znorm_drho + 0.3e1_dp*t322*t310*pnorm_drho)*t329 - &
688  0.2e1_dp*t485*t319*pnorm_drho)
689  END IF
690 
691  epsilon_crevpkzbtau = -0.2e1_dp*t263*t271
692  t799 = t264**2
693  ztau = -tau_w*t265
694  alphatau = -0.5e1_dp/0.3e1_dp*t430*ztau
695  tildeq_btau = 0.9e1_dp/0.20e2_dp*alphatau*t295 - 0.9e1_dp/ &
696  0.40e2_dp*t437*(b*alphatau*t290 + t291*alphatau)
697 
698  IF (grad_deriv >= 1 .OR. grad_deriv == -1) THEN
699  e_tau(ii) = e_tau(ii) + &
700  scale_ec*(my_rho*epsilon_crevpkzbtau*t274 + t267* &
701  (d*epsilon_crevpkzbtau*t272 - 0.3e1_dp*t268*t269/t799)) + &
702  scale_ex*t337*t490*((0.2e1_dp*t445*t303*ztau - 0.4e1_dp &
703  *t450*t452*ztau)*p + 0.292e3_dp/0.2025e4_dp*tildeq_b* &
704  tildeq_btau - 0.73e2_dp/0.4050e4_dp*tildeq_btau*t313 - &
705  0.73e2_dp/0.225e3_dp*t464*z*ztau + 0.8e1_dp/0.45e2_dp* &
706  t475*ztau)*t329
707  END IF
708  END IF
709  END IF
710  END DO
711 
712 !$OMP END DO
713 
714  END SUBROUTINE tpss_lda_calc
715 
716 END MODULE xc_tpss
bibliography
collects all references to literature in CP2K as new algorithms / method are included from literature...
Definition: bibliography.F:28
bibliography::tao2003
integer, save, public tao2003
Definition: bibliography.F:43
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
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:1040
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_tau
integer, parameter, public deriv_tau
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
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:66
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_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
xc_tpss
Calculates the tpss functional.
Definition: xc_tpss.F:17
xc_tpss::tpss_lda_eval
subroutine, public tpss_lda_eval(rho_set, deriv_set, grad_deriv, tpss_params)
evaluates the tpss functional in the spin unpolarized (lda) case
Definition: xc_tpss.F:109
xc_tpss::tpss_lda_info
subroutine, public tpss_lda_info(tpss_params, reference, shortform, needs, max_deriv)
return various information on the functional
Definition: xc_tpss.F:60