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 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 5222 : 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 2611 : CALL section_vals_val_get(tpss_params, "SCALE_C", r_val=sc)
68 2611 : CALL section_vals_val_get(tpss_params, "SCALE_X", r_val=sx)
69 :
70 2611 : IF (PRESENT(reference)) THEN
71 23 : IF (sx == 1._dp .AND. sc == 1._dp) THEN
72 23 : 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 0 : "J. Tao, J.P.Perdew, V.N.Staroverov, E.Scuseria PRL, 91, 146401 (2003)", &
76 0 : sx, sc
77 : END IF
78 : END IF
79 2611 : IF (PRESENT(shortform)) THEN
80 23 : IF (sx == 1._dp .AND. sc == 1._dp) THEN
81 23 : shortform = "TPSS meta-GGA functional (LDA)"
82 : ELSE
83 : WRITE (shortform, "(a,'sx=',f5.3,'sc=',f5.3,' (LDA)')") &
84 0 : "TPSS meta-GGA functional", &
85 0 : sx, sc
86 : END IF
87 : END IF
88 2611 : IF (PRESENT(needs)) THEN
89 2588 : needs%rho = .TRUE.
90 2588 : needs%tau = .TRUE.
91 2588 : needs%norm_drho = .TRUE.
92 : END IF
93 2611 : IF (PRESENT(max_deriv)) max_deriv = 1
94 :
95 2611 : 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 maximum degree of the derivative that should be evaluated
103 : !> \param tpss_params ...
104 : !> \author fawzi
105 : ! **************************************************************************************************
106 13544 : SUBROUTINE tpss_lda_eval(rho_set, deriv_set, grad_deriv, tpss_params)
107 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
108 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
109 : INTEGER, INTENT(in) :: grad_deriv
110 : TYPE(section_vals_type), POINTER :: tpss_params
111 :
112 : CHARACTER(len=*), PARAMETER :: routineN = 'tpss_lda_eval'
113 :
114 : INTEGER :: handle, non_coer, npoints
115 : INTEGER, DIMENSION(2, 3) :: bo
116 : REAL(kind=dp) :: epsilon_rho, epsilon_tau, scale_ec, &
117 : scale_ex
118 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
119 3386 : POINTER :: dummy, e_0, e_ndrho, e_rho, e_tau, &
120 3386 : norm_drho, rho, tau
121 : TYPE(cp_logger_type), POINTER :: logger
122 : TYPE(xc_derivative_type), POINTER :: deriv
123 :
124 3386 : CALL timeset(routineN, handle)
125 :
126 3386 : CALL cite_reference(tao2003)
127 :
128 : CALL xc_rho_set_get(rho_set, rho=rho, &
129 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
130 3386 : tau=tau, tau_cutoff=epsilon_tau)
131 3386 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
132 :
133 3386 : dummy => rho
134 :
135 3386 : e_0 => dummy
136 3386 : e_rho => dummy
137 3386 : e_ndrho => dummy
138 3386 : e_tau => dummy
139 :
140 3386 : IF (grad_deriv >= 0) THEN
141 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
142 3386 : allocate_deriv=.TRUE.)
143 3386 : CALL xc_derivative_get(deriv, deriv_data=e_0)
144 : END IF
145 3386 : IF (grad_deriv >= 1) THEN
146 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
147 3368 : allocate_deriv=.TRUE.)
148 3368 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
149 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
150 3368 : allocate_deriv=.TRUE.)
151 3368 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
152 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], &
153 3368 : allocate_deriv=.TRUE.)
154 3368 : CALL xc_derivative_get(deriv, deriv_data=e_tau)
155 : END IF
156 3386 : IF (grad_deriv > 1) THEN
157 0 : CPABORT("derivatives bigger than 1 not implemented")
158 : END IF
159 :
160 3386 : non_coer = 0
161 3386 : CALL section_vals_val_get(tpss_params, "SCALE_C", r_val=scale_ec)
162 3386 : CALL section_vals_val_get(tpss_params, "SCALE_X", r_val=scale_ex)
163 :
164 : !$OMP PARALLEL DEFAULT(NONE) &
165 : !$OMP SHARED(rho, tau, norm_drho, e_0, e_rho, e_ndrho, e_tau) &
166 : !$OMP SHARED(epsilon_rho, epsilon_tau, npoints, grad_deriv) &
167 : !$OMP SHARED(scale_ec, scale_ex) &
168 3386 : !$OMP REDUCTION(+: non_coer)
169 :
170 : CALL tpss_lda_calc(rho=rho, norm_drho=norm_drho, &
171 : tau=tau, e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, e_tau=e_tau, &
172 : grad_deriv=grad_deriv, npoints=npoints, epsilon_rho=epsilon_rho, &
173 : epsilon_tau=epsilon_tau, scale_ec=scale_ec, scale_ex=scale_ex, non_coer=non_coer)
174 :
175 : !$OMP END PARALLEL
176 :
177 3386 : logger => cp_get_default_logger()
178 : ! we could check if tau/grad were consistent, but don't do anything here
179 : IF (non_coer > 0) THEN
180 3386 : non_coer = 0
181 : END IF
182 :
183 3386 : CALL timestop(handle)
184 3386 : END SUBROUTINE tpss_lda_eval
185 :
186 : ! **************************************************************************************************
187 : !> \brief low level calculation routine for the unpolarized (lda) tpss
188 : !> \param rho ...
189 : !> \param norm_drho ...
190 : !> \param tau ...
191 : !> \param e_0 ...
192 : !> \param e_rho ...
193 : !> \param e_ndrho ...
194 : !> \param e_tau ...
195 : !> \param npoints ...
196 : !> \param grad_deriv ...
197 : !> \param epsilon_rho ...
198 : !> \param epsilon_tau ...
199 : !> \param scale_ec ...
200 : !> \param scale_ex ...
201 : !> \param non_coer ...
202 : !> \author fawzi
203 : !> \note
204 : !> maple is nice, but if you want the uman readable version of the code
205 : !> look in doc/tpss.mw
206 : ! **************************************************************************************************
207 3386 : SUBROUTINE tpss_lda_calc(rho, norm_drho, tau, e_0, e_rho, e_ndrho, e_tau, &
208 : npoints, grad_deriv, epsilon_rho, epsilon_tau, &
209 : scale_ec, scale_ex, non_coer)
210 : REAL(kind=dp), DIMENSION(*), INTENT(in) :: rho, norm_drho, tau
211 : REAL(kind=dp), DIMENSION(*), INTENT(inout) :: e_0, e_rho, e_ndrho, e_tau
212 : INTEGER, INTENT(in) :: npoints, grad_deriv
213 : REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_tau, scale_ec, &
214 : scale_ex
215 : INTEGER, INTENT(inout) :: non_coer
216 :
217 : INTEGER :: ii
218 : LOGICAL :: clipped_tau, floored_tau, t571, t639
219 : REAL(dp) :: incoming_gradient, incoming_rho, incoming_tau, partial_gradient, partial_rho, &
220 : partial_tau, tau_floor, tau_unif_inverse
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, regular_tau, rs, rs_s1, rs_s1rho, rs_s2, rs_s2rho, rsrho, t, t1, t100, t101, t111, &
230 : t12, t13, t138, t14, t140, t143, t145, t146, t147, t151, t152, t16, t161, t168, t177, &
231 : t186, t187, t189, t19, t190, t191, t193, t194, t196, t197, t198, t199, t2, t20, t201, &
232 : t202, t204, t205, t208, t209, t21, t211, t212, t213, t215, t216, t218, t219, t22, t220, &
233 : t221, t223, t224, t226, t227, t230, t231, t233, t234, t235, t238, t239, t241, t242, t243, &
234 : t245, t246, t248, t249, t252, t253, t254, t256, t26, t260, t263
235 : REAL(kind=dp) :: t264, t265, t267, t268, t269, t27, t271, t272, t274, t275, t276, t277, &
236 : t278, t279, t28, t280, t281, t284, t286, t29, t290, t291, t293, t294, t295, t299, t3, &
237 : t301, t302, t303, t305, t307, t310, t313, t316, t319, t322, t325, t327, t328, t329, t331, &
238 : t337, t340, t343, t344, t35, t351, t36, t370, t371, t376, t383, t385, t386, t39, t390, &
239 : t391, t395, t396, t398, t4, t403, t404, t406, t41, t410, t411, t419, t42, t437, t445, &
240 : t450, t452, t464, t472, t475, t485, t489, t49, t490, t5, t505, t513, t517, t536, t541, &
241 : t542, 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 3386 : kappa = 0.804e0_dp
249 3386 : beta = 0.66725e-1_dp
250 3386 : mu = 0.21951e0_dp
251 3386 : gamma_var = (0.1e1_dp - LOG(0.2e1_dp))/pi**2
252 3386 : b = 0.4e0_dp
253 3386 : c = 0.159096e1_dp
254 3386 : e_var = 0.1537e1_dp
255 3386 : d = 0.28e1_dp
256 3386 : p_1 = 0.10e1_dp
257 3386 : A_1 = 0.31091e-1_dp
258 3386 : alpha_1_1 = 0.21370e0_dp
259 3386 : beta_1_1 = 0.75957e1_dp
260 3386 : beta_2_1 = 0.35876e1_dp
261 3386 : beta_3_1 = 0.16382e1_dp
262 3386 : beta_4_1 = 0.49294e0_dp
263 3386 : p_2 = 0.10e1_dp
264 3386 : A_2 = 0.15545e-1_dp
265 3386 : alpha_1_2 = 0.20548e0_dp
266 3386 : beta_1_2 = 0.141189e2_dp
267 3386 : beta_2_2 = 0.61977e1_dp
268 3386 : beta_3_2 = 0.33662e1_dp
269 3386 : beta_4_2 = 0.62517e0_dp
270 3386 : p_3 = 0.10e1_dp
271 :
272 3386 : tau_floor = MAX(epsilon_tau, EPSILON(0.0_dp))
273 3386 : t1 = 3._dp**(0.1e1_dp/0.3e1_dp)
274 3386 : t2 = 4._dp**(0.1e1_dp/0.3e1_dp)
275 3386 : t3 = t2**2
276 3386 : t4 = t1*t3
277 3386 : t5 = 2._dp**(0.1e1_dp/0.3e1_dp)
278 3386 : t6 = 0.1e1_dp/pi
279 3386 : t12 = t5**2
280 :
281 3386 : !$OMP DO
282 :
283 : DO ii = 1, npoints
284 : ! Fourier interpolation can give negative tau even at finite density.
285 : ! Continue the energy through tau=0; discarding these points makes it
286 : ! discontinuous and invalidates the potential and geometry derivatives.
287 126007345 : floored_tau = tau(ii) <= tau_floor
288 126007345 : regular_tau = MAX(tau(ii), tau_floor)
289 126007345 : my_tau = regular_tau
290 126007345 : my_rho = rho(ii)
291 126007345 : IF (my_rho > epsilon_rho) THEN
292 121554717 : my_ndrho = norm_drho(ii)
293 121554717 : IF (grad_deriv > 0) THEN
294 119403449 : incoming_rho = e_rho(ii)
295 119403449 : incoming_gradient = e_ndrho(ii)
296 119403449 : incoming_tau = e_tau(ii)
297 : END IF
298 :
299 121554717 : t7 = 0.1e1_dp/my_rho
300 121554717 : t254 = my_ndrho**2
301 121554717 : tau_w = t254*t7/0.8e1_dp
302 :
303 121554717 : clipped_tau = my_tau < tau_w
304 121554717 : IF (clipped_tau) THEN
305 : ! Project the gradient and kinetic density onto z=1.
306 2724599 : m = 0.5_dp*t254 + 4.0_dp*my_rho*my_tau
307 2724599 : my_tau = m/8._dp/my_rho
308 2724599 : my_ndrho = SQRT(m)
309 2724599 : t254 = m
310 2724599 : tau_w = my_tau
311 2724599 : non_coer = non_coer + 1
312 : END IF
313 :
314 121554717 : t9 = (t6*t7)**(0.1e1_dp/0.3e1_dp)
315 121554717 : rs_s1 = t4*t5*t9/0.4e1_dp
316 121554717 : phi_s1 = t12/0.2e1_dp
317 121554717 : t13 = t1*t12
318 121554717 : t14 = pi**2
319 121554717 : t16 = (t14*my_rho)**(0.1e1_dp/0.3e1_dp)
320 121554717 : k_f_s1 = t13*t16/0.2e1_dp
321 121554717 : t19 = SQRT(k_f_s1*t6)
322 121554717 : k_s_s1 = 0.2e1_dp*t19
323 121554717 : t20 = 0.1e1_dp/phi_s1
324 121554717 : t21 = my_ndrho*t20
325 121554717 : t22 = 0.1e1_dp/k_s_s1
326 121554717 : t_s1 = t21*t22*t7/0.2e1_dp
327 121554717 : rs_s2 = rs_s1
328 121554717 : phi_s2 = phi_s1
329 121554717 : t26 = SQRT(k_f_s1*t6)
330 121554717 : k_s_s2 = 0.2e1_dp*t26
331 121554717 : t27 = 0.1e1_dp/phi_s2
332 121554717 : t28 = my_ndrho*t27
333 121554717 : t29 = 0.1e1_dp/k_s_s2
334 121554717 : t_s2 = t28*t29*t7/0.2e1_dp
335 121554717 : t35 = 0.1e1_dp/A_1
336 121554717 : t36 = SQRT(rs_s2)
337 121554717 : t39 = t36*rs_s2
338 121554717 : t41 = p_1 + 0.1e1_dp
339 121554717 : t42 = rs_s2**t41
340 : t49 = LOG(0.1e1_dp + t35/(beta_1_1*t36 + beta_2_1*rs_s2 + &
341 : beta_3_1*t39 + beta_4_1*t42)/0.2e1_dp)
342 121554717 : t55 = SQRT(rs_s1)
343 121554717 : t58 = t55*rs_s1
344 121554717 : t60 = rs_s1**t41
345 : t67 = LOG(0.1e1_dp + t35/(beta_1_1*t55 + beta_2_1*rs_s1 + &
346 : beta_3_1*t58 + beta_4_1*t60)/0.2e1_dp)
347 121554717 : t71 = 0.1e1_dp + alpha_1_2*rs_s2
348 121554717 : t73 = 0.1e1_dp/A_2
349 121554717 : t77 = p_2 + 0.1e1_dp
350 121554717 : t78 = rs_s2**t77
351 121554717 : t79 = beta_4_2*t78
352 121554717 : t80 = beta_1_2*t36 + beta_2_2*rs_s2 + beta_3_2*t39 + t79
353 121554717 : t84 = 0.1e1_dp + t73/t80/0.2e1_dp
354 121554717 : t85 = LOG(t84)
355 121554717 : e_c_u_1_s2 = -0.2e1_dp*A_2*t71*t85
356 121554717 : t89 = 0.1e1_dp + alpha_1_2*rs_s1
357 121554717 : t94 = rs_s1**t77
358 121554717 : t95 = beta_4_2*t94
359 121554717 : t96 = beta_1_2*t55 + beta_2_2*rs_s1 + beta_3_2*t58 + t95
360 121554717 : t100 = 0.1e1_dp + t73/t96/0.2e1_dp
361 121554717 : t101 = LOG(t100)
362 121554717 : e_c_u_1_s1 = -0.2e1_dp*A_2*t89*t101
363 121554717 : t111 = p_3 + 1._dp
364 121554717 : rs = t4*t9/0.4e1_dp
365 121554717 : t138 = 0.1e1_dp + alpha_1_1*rs
366 121554717 : t140 = SQRT(rs)
367 121554717 : t143 = t140*rs
368 121554717 : t145 = rs**t41
369 121554717 : t146 = beta_4_1*t145
370 121554717 : t147 = beta_1_1*t140 + beta_2_1*rs + beta_3_1*t143 + t146
371 121554717 : t151 = 0.1e1_dp + t35/t147/0.2e1_dp
372 121554717 : t152 = LOG(t151)
373 121554717 : e_c_u_0 = -0.2e1_dp*A_1*t138*t152
374 121554717 : t161 = rs**t77
375 : t168 = LOG(0.1e1_dp + t73/(beta_1_2*t140 + beta_2_2*rs + &
376 : beta_3_2*t143 + beta_4_2*t161)/0.2e1_dp)
377 121554717 : t177 = rs**t111
378 121554717 : t186 = 0.1e1_dp/gamma_var
379 121554717 : t187 = beta*t186
380 121554717 : t189 = phi_s1**2
381 121554717 : t190 = t189*phi_s1
382 121554717 : t191 = 0.1e1_dp/t190
383 121554717 : t193 = EXP(-e_c_u_1_s1*t186*t191)
384 121554717 : t194 = t193 - 0.1e1_dp
385 121554717 : A_s1 = t187/t194
386 121554717 : t196 = gamma_var*t190
387 121554717 : t197 = t_s1**2
388 121554717 : t198 = A_s1*t197
389 121554717 : t199 = 0.1e1_dp + t198
390 121554717 : t201 = A_s1**2
391 121554717 : t202 = t197**2
392 121554717 : t204 = 0.1e1_dp + t198 + t201*t202
393 121554717 : t205 = 0.1e1_dp/t204
394 121554717 : t208 = 0.1e1_dp + t187*t197*t199*t205
395 121554717 : t209 = LOG(t208)
396 121554717 : epsilon_cGGA_1_0 = e_c_u_1_s1 + t196*t209
397 121554717 : t211 = phi_s2**2
398 121554717 : t212 = t211*phi_s2
399 121554717 : t213 = 0.1e1_dp/t212
400 121554717 : t215 = EXP(-e_c_u_1_s2*t186*t213)
401 121554717 : t216 = t215 - 0.1e1_dp
402 121554717 : A_s2 = t187/t216
403 121554717 : t218 = gamma_var*t212
404 121554717 : t219 = t_s2**2
405 121554717 : t220 = A_s2*t219
406 121554717 : t221 = t220 + 0.1e1_dp
407 121554717 : t223 = A_s2**2
408 121554717 : t224 = t219**2
409 121554717 : t226 = 0.1e1_dp + t220 + t223*t224
410 121554717 : t227 = 0.1e1_dp/t226
411 121554717 : t230 = 0.1e1_dp + t187*t219*t221*t227
412 121554717 : t231 = LOG(t230)
413 121554717 : epsilon_cGGA_0_1 = e_c_u_1_s2 + t218*t231
414 121554717 : t233 = SQRT(t1*t16*t6)
415 121554717 : k_s = 0.2e1_dp*t233
416 121554717 : t234 = 0.1e1_dp/k_s
417 121554717 : t235 = my_ndrho*t234
418 121554717 : t = t235*t7/0.2e1_dp
419 121554717 : t238 = EXP(-e_c_u_0*t186)
420 121554717 : t239 = -0.1e1_dp + t238
421 121554717 : A = t187/t239
422 121554717 : t241 = t**2
423 121554717 : t242 = A*t241
424 121554717 : t243 = 0.1e1_dp + t242
425 121554717 : t245 = A**2
426 121554717 : t246 = t241**2
427 121554717 : t248 = 0.1e1_dp + t242 + t245*t246
428 121554717 : t249 = 0.1e1_dp/t248
429 121554717 : t252 = 0.1e1_dp + t187*t241*t243*t249
430 121554717 : t253 = LOG(t252)
431 121554717 : epsilon_cGGA = e_c_u_0 + gamma_var*t253
432 121554717 : ma = MAX(epsilon_cGGA_1_0, epsilon_cGGA)
433 121554717 : mb = MAX(epsilon_cGGA_0_1, epsilon_cGGA)
434 121554717 : t256 = tau_w**2
435 121554717 : t260 = ma/0.2e1_dp + mb/0.2e1_dp
436 121554717 : t263 = 0.53e0_dp*epsilon_cGGA*t256 - 0.153e1_dp*t256*t260
437 121554717 : t264 = my_tau**2
438 121554717 : t265 = 0.1e1_dp/t264
439 121554717 : epsilon_cRevPKZB = epsilon_cGGA + t263*t265
440 121554717 : t267 = my_rho*epsilon_cRevPKZB
441 121554717 : t268 = d*epsilon_cRevPKZB
442 121554717 : t269 = t256*tau_w
443 121554717 : t271 = 0.1e1_dp/t264/my_tau
444 121554717 : t272 = t269*t271
445 121554717 : t274 = 0.1e1_dp + t268*t272
446 121554717 : t275 = t254*t1
447 121554717 : t276 = t14**(0.1e1_dp/0.3e1_dp)
448 121554717 : t277 = t276**2
449 121554717 : t278 = 0.1e1_dp/t277
450 121554717 : t279 = my_rho**2
451 121554717 : t280 = my_rho**(0.1e1_dp/0.3e1_dp)
452 121554717 : t281 = t280**2
453 121554717 : t284 = t278/t281/t279
454 121554717 : p = t275*t284/0.12e2_dp
455 121554717 : t286 = 0.1e1_dp/my_tau
456 121554717 : z = tau_w*t286
457 : ! Use the finite stationary-density limit of (5/3)*p*(1/z-1).
458 121554717 : tau_unif_inverse = 10.0_dp/9.0_dp*t1*t278/(t281*my_rho)
459 121554717 : alpha = (my_tau - tau_w)*tau_unif_inverse
460 121554717 : t290 = alpha - 0.1e1_dp
461 121554717 : t291 = b*alpha
462 121554717 : t293 = 0.1e1_dp + t291*t290
463 121554717 : t294 = SQRT(t293)
464 121554717 : t295 = 0.1e1_dp/t294
465 121554717 : tildeq_b = 0.9e1_dp/0.20e2_dp*t290*t295 + 0.2e1_dp/0.3e1_dp*p
466 121554717 : t299 = z**2
467 121554717 : t301 = 0.1e1_dp + t299
468 121554717 : t302 = t301**2
469 121554717 : t303 = 0.1e1_dp/t302
470 121554717 : t305 = 0.10e2_dp/0.81e2_dp + c*t299*t303
471 121554717 : t307 = tildeq_b**2
472 121554717 : t310 = p**2
473 121554717 : t313 = SQRT(0.18e2_dp*t299 + 0.50e2_dp*t310)
474 121554717 : t316 = 0.1e1_dp/kappa
475 121554717 : t319 = SQRT(e_var)
476 121554717 : t322 = e_var*mu
477 : t325 = t305*p + 0.146e3_dp/0.2025e4_dp*t307 - 0.73e2_dp/ &
478 : 0.4050e4_dp*tildeq_b*t313 + 0.100e3_dp/0.6561e4_dp*t316* &
479 121554717 : t310 + 0.4e1_dp/0.45e2_dp*t319*t299 + t322*t310*p
480 121554717 : t327 = 0.1e1_dp + t319*p
481 121554717 : t328 = t327**2
482 121554717 : t329 = 0.1e1_dp/t328
483 121554717 : t331 = 0.1e1_dp + t325*t329*t316
484 121554717 : Fx = 0.1e1_dp + kappa - kappa/t331
485 121554717 : ex_unif = -0.3e1_dp/0.4e1_dp*t1*t16*t6
486 121554717 : t337 = my_rho*ex_unif
487 :
488 : e_0(ii) = e_0(ii) + &
489 121554717 : scale_ec*t267*t274 + scale_ex*t337*Fx
490 :
491 121554717 : IF (grad_deriv > 0) THEN
492 119403449 : t340 = t9**2
493 119403449 : t343 = 0.1e1_dp/t279
494 119403449 : t344 = 0.1e1_dp/t340*t6*t343
495 119403449 : rsrho = -t4*t344/0.12e2_dp
496 119403449 : t351 = t147**2
497 : e_c_u_0rho = -0.2e1_dp*A_1*alpha_1_1*rsrho*t152 + t138/ &
498 : t351*(beta_1_1/t140*rsrho/0.2e1_dp + beta_2_1*rsrho + &
499 : 0.3e1_dp/0.2e1_dp*beta_3_1*t140*rsrho + t146*t41*rsrho/ &
500 119403449 : rs)/t151
501 119403449 : t370 = t16**2
502 119403449 : t371 = 0.1e1_dp/t370
503 119403449 : t376 = k_s**2
504 : trho = -my_ndrho/t376*t7/t233*t1*t371*t14*t6 &
505 119403449 : /0.6e1_dp - t235*t343/0.2e1_dp
506 119403449 : t383 = gamma_var**2
507 119403449 : t385 = beta/t383
508 119403449 : t386 = t239**2
509 119403449 : Arho = t385/t386*e_c_u_0rho*t238
510 119403449 : t390 = t187*t
511 119403449 : t391 = t243*t249
512 119403449 : t395 = Arho*t241
513 119403449 : t396 = A*t
514 119403449 : t398 = 0.2e1_dp*t396*trho
515 119403449 : t403 = t187*t241
516 119403449 : t404 = t248**2
517 119403449 : t406 = t243/t404
518 119403449 : t410 = t241*t
519 119403449 : t411 = t245*t410
520 119403449 : t419 = 0.1e1_dp/t252
521 : epsilon_cGGArho = e_c_u_0rho + gamma_var*(0.2e1_dp*t390*t391 &
522 : *trho + t187*t241*(t395 + t398)*t249 - t403*t406*(t395 + &
523 119403449 : t398 + 0.2e1_dp*A*t246*Arho + 0.4e1_dp*t411*trho))*t419
524 119403449 : tau_wrho = -t254*t343/0.8e1_dp
525 119403449 : prho = -0.2e1_dp/0.9e1_dp*t275*t278/t281/t279/my_rho
526 119403449 : zrho = tau_wrho*t286
527 119403449 : alpharho = 5.0_dp/3.0_dp*(p - alpha)*t7
528 119403449 : t437 = t290/t294/t293
529 : tildeq_brho = 0.9e1_dp/0.20e2_dp*alpharho*t295 - 0.9e1_dp/ &
530 : 0.40e2_dp*t437*(b*alpharho*t290 + t291*alpharho) + &
531 119403449 : 0.2e1_dp/0.3e1_dp*prho
532 119403449 : t445 = c*z
533 119403449 : t450 = c*t299*z
534 119403449 : t452 = 0.1e1_dp/t302/t301
535 : ! The associated first-gradient factors vanish at t313=0.
536 119403449 : t464 = 0.0_dp
537 119403449 : IF (t313 > 0.0_dp) t464 = tildeq_b/t313
538 119403449 : t472 = t316*p
539 119403449 : t475 = t319*z
540 119403449 : t485 = t325/t328/t327
541 119403449 : t489 = t331**2
542 119403449 : t490 = 0.1e1_dp/t489
543 119403449 : rs_s1rho = -t4*t5*t344/0.12e2_dp
544 119403449 : k_f_s1rho = t13*t371*t14/0.6e1_dp
545 119403449 : t505 = k_s_s1**2
546 : t_s1rho = -t21/t505*t7/t19*k_f_s1rho*t6/0.2e1_dp - t21 &
547 119403449 : *t22*t343/0.2e1_dp
548 119403449 : t513 = A_2*alpha_1_2
549 119403449 : t517 = t96**2
550 : e_c_u_1_s1rho = -0.2e1_dp*t513*rs_s1rho*t101 + t89/t517* &
551 : (beta_1_2/t55*rs_s1rho/0.2e1_dp + beta_2_2*rs_s1rho + &
552 : 0.3e1_dp/0.2e1_dp*beta_3_2*t55*rs_s1rho + t95*t77* &
553 119403449 : rs_s1rho/rs_s1)/t100
554 119403449 : t536 = t194**2
555 119403449 : A_s1rho = t385/t536*e_c_u_1_s1rho*t191*t193
556 119403449 : t541 = t187*t_s1
557 119403449 : t542 = t199*t205
558 119403449 : t546 = A_s1rho*t197
559 119403449 : t547 = A_s1*t_s1
560 119403449 : t549 = 0.2e1_dp*t547*t_s1rho
561 119403449 : t554 = t187*t197
562 119403449 : t555 = t204**2
563 119403449 : t557 = t199/t555
564 119403449 : t561 = t197*t_s1
565 119403449 : t562 = t201*t561
566 119403449 : t569 = 0.1e1_dp/t208
567 119403449 : t571 = epsilon_cGGA < epsilon_cGGA_1_0
568 119403449 : IF (t571) THEN
569 : marho = e_c_u_1_s1rho + t196*(0.2e1_dp*t541*t542 &
570 : *t_s1rho + t187*t197*(t546 + t549)*t205 - t554*t557*(t546 &
571 : + t549 + 0.2e1_dp*A_s1*t202*A_s1rho + 0.4e1_dp*t562* &
572 119330471 : t_s1rho))*t569
573 : ELSE
574 : marho = epsilon_cGGArho
575 : END IF
576 119403449 : rs_s2rho = rs_s1rho
577 119403449 : t574 = k_s_s2**2
578 : t_s2rho = -t28/t574*t7/t26*k_f_s1rho*t6/0.2e1_dp - t28 &
579 119403449 : *t29*t343/0.2e1_dp
580 119403449 : t585 = t80**2
581 : e_c_u_1_s2rho = -0.2e1_dp*t513*rs_s2rho*t85 + t71/t585*( &
582 : beta_1_2/t36*rs_s2rho/0.2e1_dp + beta_2_2*rs_s2rho + &
583 : 0.3e1_dp/0.2e1_dp*beta_3_2*t36*rs_s2rho + t79*t77* &
584 119403449 : rs_s2rho/rs_s2)/t84
585 119403449 : t604 = t216**2
586 119403449 : A_s2rho = t385/t604*e_c_u_1_s2rho*t213*t215
587 119403449 : t609 = t187*t_s2
588 119403449 : t610 = t221*t227
589 119403449 : t614 = A_s2rho*t219
590 119403449 : t615 = A_s2*t_s2
591 119403449 : t617 = 0.2e1_dp*t615*t_s2rho
592 119403449 : t622 = t187*t219
593 119403449 : t623 = t226**2
594 119403449 : t625 = t221/t623
595 119403449 : t629 = t219*t_s2
596 119403449 : t630 = t223*t629
597 119403449 : t637 = 0.1e1_dp/t230
598 119403449 : t639 = epsilon_cGGA < epsilon_cGGA_0_1
599 : IF (t639) THEN
600 : mbrho = e_c_u_1_s2rho + t218*(0.2e1_dp*t609*t610 &
601 : *t_s2rho + t187*t219*(t614 + t617)*t227 - t622*t625*(t614 &
602 : + t617 + 0.2e1_dp*A_s2*t224*A_s2rho + 0.4e1_dp*t630* &
603 119403449 : t_s2rho))*t637
604 : ELSE
605 : mbrho = epsilon_cGGArho
606 : END IF
607 119403449 : t642 = epsilon_cGGA*tau_w
608 119403449 : t645 = tau_w*t260
609 : epsilon_cRevPKZBrho = epsilon_cGGArho + (0.53e0_dp* &
610 : epsilon_cGGArho*t256 + 0.106e1_dp*t642*tau_wrho - 0.306e1_dp* &
611 : t645*tau_wrho - 0.153e1_dp*t256*(marho/0.2e1_dp + mbrho/ &
612 119403449 : 0.2e1_dp))*t265
613 119403449 : t659 = t256*t271
614 :
615 : e_rho(ii) = e_rho(ii) + &
616 : scale_ec*(epsilon_cRevPKZB*t274 + my_rho* &
617 : epsilon_cRevPKZBrho*t274 + t267*(d*epsilon_cRevPKZBrho*t272 &
618 : + 0.3e1_dp*t268*t659*tau_wrho)) + scale_ex*(ex_unif*Fx - &
619 : my_rho*pi*t1*t371*Fx/0.4e1_dp + t337* &
620 : t490*(((0.2e1_dp*t445*t303*zrho - 0.4e1_dp*t450*t452* &
621 : zrho)*p + t305*prho + 0.292e3_dp/0.2025e4_dp*tildeq_b* &
622 : tildeq_brho - 0.73e2_dp/0.4050e4_dp*tildeq_brho*t313 - &
623 : 0.73e2_dp/0.8100e4_dp*t464*(0.36e2_dp*z*zrho + 0.100e3_dp &
624 : *p*prho) + 0.200e3_dp/0.6561e4_dp*t472*prho + 0.8e1_dp/ &
625 : 0.45e2_dp*t475*zrho + 0.3e1_dp*t322*t310*prho)*t329 - &
626 119403449 : 0.2e1_dp*t485*t319*prho))
627 :
628 119403449 : tnorm_drho = t234*t7/0.2e1_dp
629 : Hnorm_drho = gamma_var*(0.2e1_dp*t390*t391*tnorm_drho + &
630 : 0.2e1_dp*t187*t410*A*tnorm_drho*t249 - t403*t406*( &
631 119403449 : 0.2e1_dp*t396*tnorm_drho + 0.4e1_dp*t411*tnorm_drho))*t419
632 119403449 : tau_wnorm_drho = my_ndrho*t7/0.4e1_dp
633 119403449 : pnorm_drho = my_ndrho*t1*t284/0.6e1_dp
634 119403449 : znorm_drho = tau_wnorm_drho*t286
635 119403449 : alphanorm_drho = -5.0_dp/3.0_dp*pnorm_drho
636 : tildeq_bnorm_drho = 0.9e1_dp/0.20e2_dp*alphanorm_drho*t295 - &
637 : 0.9e1_dp/0.40e2_dp*t437*(b*alphanorm_drho*t290 + t291* &
638 119403449 : alphanorm_drho) + 0.2e1_dp/0.3e1_dp*pnorm_drho
639 119403449 : t_s1norm_drho = t20*t22*t7/0.2e1_dp
640 119403449 : IF (t571) THEN
641 : manorm_drho = t196*(0.2e1_dp*t541*t542* &
642 : t_s1norm_drho + 0.2e1_dp*t187*t561*A_s1*t_s1norm_drho*t205 &
643 : - t554*t557*(0.2e1_dp*t547*t_s1norm_drho + 0.4e1_dp*t562 &
644 119330471 : *t_s1norm_drho))*t569
645 : ELSE
646 : manorm_drho = Hnorm_drho
647 : END IF
648 119403449 : t_s2norm_drho = t27*t29*t7/0.2e1_dp
649 : IF (t639) THEN
650 : mbnorm_drho = t218*(0.2e1_dp*t609*t610* &
651 : t_s2norm_drho + 0.2e1_dp*t187*t629*A_s2*t_s2norm_drho*t227 &
652 : - t622*t625*(0.2e1_dp*t615*t_s2norm_drho + 0.4e1_dp*t630 &
653 119403449 : *t_s2norm_drho))*t637
654 : ELSE
655 : mbnorm_drho = Hnorm_drho
656 : END IF
657 : epsilon_cRevPKZBnorm_drho = Hnorm_drho + (0.53e0_dp*Hnorm_drho* &
658 : t256 + 0.106e1_dp*t642*tau_wnorm_drho - 0.306e1_dp*t645* &
659 : tau_wnorm_drho - 0.153e1_dp*t256*(manorm_drho/0.2e1_dp + &
660 119403449 : mbnorm_drho/0.2e1_dp))*t265
661 :
662 : e_ndrho(ii) = e_ndrho(ii) + &
663 : scale_ec*(my_rho*epsilon_cRevPKZBnorm_drho* &
664 : t274 + t267*(d*epsilon_cRevPKZBnorm_drho*t272 + 0.3e1_dp* &
665 : t268*t659*tau_wnorm_drho)) + scale_ex*t337*t490*((( &
666 : 0.2e1_dp*t445*t303*znorm_drho - 0.4e1_dp*t450*t452* &
667 : znorm_drho)*p + t305*pnorm_drho + 0.292e3_dp/0.2025e4_dp* &
668 : tildeq_b*tildeq_bnorm_drho - 0.73e2_dp/0.4050e4_dp* &
669 : tildeq_bnorm_drho*t313 - 0.73e2_dp/0.8100e4_dp*t464*( &
670 : 0.36e2_dp*z*znorm_drho + 0.100e3_dp*p*pnorm_drho) + &
671 : 0.200e3_dp/0.6561e4_dp*t472*pnorm_drho + 0.8e1_dp/0.45e2_dp &
672 : *t475*znorm_drho + 0.3e1_dp*t322*t310*pnorm_drho)*t329 - &
673 119403449 : 0.2e1_dp*t485*t319*pnorm_drho)
674 :
675 119403449 : epsilon_cRevPKZBtau = -0.2e1_dp*t263*t271
676 119403449 : t799 = t264**2
677 119403449 : ztau = -tau_w*t265
678 119403449 : alphatau = tau_unif_inverse
679 : tildeq_btau = 0.9e1_dp/0.20e2_dp*alphatau*t295 - 0.9e1_dp/ &
680 119403449 : 0.40e2_dp*t437*(b*alphatau*t290 + t291*alphatau)
681 :
682 : e_tau(ii) = e_tau(ii) + &
683 : scale_ec*(my_rho*epsilon_cRevPKZBtau*t274 + t267* &
684 : (d*epsilon_cRevPKZBtau*t272 - 0.3e1_dp*t268*t269/t799)) + &
685 : scale_ex*t337*t490*((0.2e1_dp*t445*t303*ztau - 0.4e1_dp &
686 : *t450*t452*ztau)*p + 0.292e3_dp/0.2025e4_dp*tildeq_b* &
687 : tildeq_btau - 0.73e2_dp/0.4050e4_dp*tildeq_btau*t313 - &
688 : 0.73e2_dp/0.225e3_dp*t464*z*ztau + 0.8e1_dp/0.45e2_dp* &
689 119403449 : t475*ztau)*t329
690 119403449 : IF (clipped_tau) THEN
691 : ! Differentiate the same projected gradient and kinetic density
692 : ! used by the energy, preserving other functional contributions.
693 2712252 : partial_rho = e_rho(ii) - incoming_rho
694 2712252 : partial_gradient = e_ndrho(ii) - incoming_gradient
695 2712252 : partial_tau = e_tau(ii) - incoming_tau
696 : e_rho(ii) = incoming_rho + partial_rho + &
697 : 2.0_dp*regular_tau/my_ndrho*partial_gradient - &
698 2712252 : norm_drho(ii)**2/(16.0_dp*my_rho**2)*partial_tau
699 : e_ndrho(ii) = incoming_gradient + &
700 : 0.5_dp*norm_drho(ii)/my_ndrho*partial_gradient + &
701 2712252 : norm_drho(ii)/(8.0_dp*my_rho)*partial_tau
702 2712252 : e_tau(ii) = incoming_tau + 2.0_dp*my_rho/my_ndrho*partial_gradient + 0.5_dp*partial_tau
703 : END IF
704 119403449 : IF (floored_tau) e_tau(ii) = incoming_tau
705 : END IF
706 : END IF
707 : END DO
708 :
709 : !$OMP END DO
710 :
711 3386 : END SUBROUTINE tpss_lda_calc
712 :
713 : END MODULE xc_tpss
|
