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 Calculate the Perdew Correlation from 1986
10 : !> \par History
11 : !> JGH (26.02.2003) : OpenMP enabled
12 : !> fawzi (04.2004) : adapted to the new xc interface
13 : !> \author JGH (03.03.2002)
14 : ! **************************************************************************************************
15 : MODULE xc_perdew86
16 :
17 : USE input_section_types, ONLY: section_vals_type
18 : USE kinds, ONLY: dp
19 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
20 : deriv_rho
21 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
22 : xc_dset_get_derivative
23 : USE xc_derivative_types, ONLY: xc_derivative_get,&
24 : xc_derivative_type
25 : USE xc_functionals_utilities, ONLY: calc_rs_pw,&
26 : set_util
27 : USE xc_input_constants, ONLY: pz_orig
28 : USE xc_perdew_zunger, ONLY: pz_lda_eval
29 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
30 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
31 : xc_rho_set_type
32 : #include "../base/base_uses.f90"
33 :
34 : IMPLICIT NONE
35 :
36 : PRIVATE
37 :
38 : ! *** Global parameters ***
39 : REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
40 : f23 = 2.0_dp*f13, &
41 : f43 = 4.0_dp*f13, &
42 : f53 = 5.0_dp*f13, &
43 : f76 = 7.0_dp/6.0_dp, &
44 : frs = 1.6119919540164696407_dp, &
45 : fpe = 0.19199566167376364_dp
46 :
47 : PUBLIC :: p86_lda_info, p86_lda_eval
48 :
49 : REAL(KIND=dp) :: eps_rho
50 : LOGICAL :: debug_flag
51 :
52 : REAL(KIND=dp), PARAMETER :: a = 0.023266_dp, &
53 : b = 7.389e-6_dp, &
54 : c = 8.723_dp, &
55 : d = 0.472_dp, &
56 : pc1 = 0.001667_dp, &
57 : pc2 = 0.002568_dp, &
58 : pci = pc1 + pc2
59 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_perdew86'
60 :
61 : CONTAINS
62 :
63 : ! **************************************************************************************************
64 : !> \brief ...
65 : !> \param cutoff ...
66 : !> \param debug ...
67 : ! **************************************************************************************************
68 180 : SUBROUTINE p86_init(cutoff, debug)
69 :
70 : REAL(KIND=dp), INTENT(IN) :: cutoff
71 : LOGICAL, INTENT(IN), OPTIONAL :: debug
72 :
73 180 : eps_rho = cutoff
74 180 : CALL set_util(cutoff)
75 :
76 180 : IF (PRESENT(debug)) THEN
77 0 : debug_flag = debug
78 : ELSE
79 180 : debug_flag = .FALSE.
80 : END IF
81 :
82 180 : END SUBROUTINE p86_init
83 :
84 : ! **************************************************************************************************
85 : !> \brief ...
86 : !> \param reference ...
87 : !> \param shortform ...
88 : !> \param needs ...
89 : !> \param max_deriv ...
90 : ! **************************************************************************************************
91 215 : SUBROUTINE p86_lda_info(reference, shortform, needs, max_deriv)
92 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
93 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
94 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
95 :
96 215 : IF (PRESENT(reference)) THEN
97 7 : reference = "J. P. Perdew, Phys. Rev. B, 33, 8822 (1986) {LDA version}"
98 : END IF
99 215 : IF (PRESENT(shortform)) THEN
100 7 : shortform = "Perdew 1986 correlation energy functional {LDA}"
101 : END IF
102 215 : IF (PRESENT(needs)) THEN
103 208 : needs%rho = .TRUE.
104 208 : needs%norm_drho = .TRUE.
105 : END IF
106 215 : IF (PRESENT(max_deriv)) max_deriv = 3
107 :
108 215 : END SUBROUTINE p86_lda_info
109 :
110 : ! **************************************************************************************************
111 : !> \brief ...
112 : !> \param rho_set ...
113 : !> \param deriv_set ...
114 : !> \param order ...
115 : !> \param p86_params ...
116 : ! **************************************************************************************************
117 180 : SUBROUTINE p86_lda_eval(rho_set, deriv_set, order, p86_params)
118 :
119 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
120 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
121 : INTEGER, INTENT(IN) :: order
122 : TYPE(section_vals_type), POINTER :: p86_params
123 :
124 : CHARACTER(len=*), PARAMETER :: routineN = 'p86_lda_eval'
125 :
126 : INTEGER :: handle, m, npoints
127 : INTEGER, DIMENSION(2, 3) :: bo
128 : REAL(KIND=dp) :: drho_cutoff, rho_cutoff
129 180 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rs
130 180 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
131 180 : e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
132 180 : e_rho_rho_rho, grho, rho
133 : TYPE(xc_derivative_type), POINTER :: deriv
134 :
135 180 : CALL timeset(routineN, handle)
136 180 : NULLIFY (rho, e_0, e_rho, e_ndrho, &
137 180 : e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
138 180 : e_rho_rho_rho, e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho)
139 :
140 : ! calculate the perdew_zunger correlation
141 180 : CALL pz_lda_eval(pz_orig, rho_set, deriv_set, order, p86_params)
142 :
143 : CALL xc_rho_set_get(rho_set, rho=rho, &
144 : norm_drho=grho, local_bounds=bo, rho_cutoff=rho_cutoff, &
145 180 : drho_cutoff=drho_cutoff)
146 180 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
147 180 : CALL p86_init(rho_cutoff)
148 180 : m = ABS(order)
149 :
150 540 : ALLOCATE (rs(npoints))
151 :
152 180 : CALL calc_rs_pw(rho, rs, npoints)
153 180 : IF (order >= 0) THEN
154 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
155 180 : allocate_deriv=.TRUE.)
156 180 : CALL xc_derivative_get(deriv, deriv_data=e_0)
157 :
158 180 : CALL p86_u_0(rho, rs, grho, e_0, npoints)
159 : END IF
160 180 : IF (order >= 1 .OR. order == -1) THEN
161 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
162 180 : allocate_deriv=.TRUE.)
163 180 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
164 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
165 180 : allocate_deriv=.TRUE.)
166 180 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
167 :
168 : CALL p86_u_1(rho, grho, rs, e_rho, &
169 180 : e_ndrho, npoints)
170 : END IF
171 180 : IF (order >= 2 .OR. order == -2) THEN
172 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
173 0 : allocate_deriv=.TRUE.)
174 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
175 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
176 0 : allocate_deriv=.TRUE.)
177 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
178 : deriv => xc_dset_get_derivative(deriv_set, &
179 0 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
180 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
181 :
182 : CALL p86_u_2(rho, grho, rs, e_rho_rho, &
183 0 : e_rho_ndrho, e_ndrho_ndrho, npoints)
184 : END IF
185 180 : IF (order >= 3 .OR. order == -3) THEN
186 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
187 0 : allocate_deriv=.TRUE.)
188 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
189 : deriv => xc_dset_get_derivative(deriv_set, &
190 0 : [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.TRUE.)
191 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
192 : deriv => xc_dset_get_derivative(deriv_set, &
193 0 : [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
194 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
195 : deriv => xc_dset_get_derivative(deriv_set, &
196 0 : [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
197 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
198 :
199 : CALL p86_u_3(rho, grho, rs, e_rho_rho_rho, &
200 : e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
201 0 : npoints)
202 : END IF
203 180 : IF (order > 3 .OR. order < -3) THEN
204 0 : CPABORT("derivatives bigger than 3 not implemented")
205 : END IF
206 180 : DEALLOCATE (rs)
207 180 : CALL timestop(handle)
208 :
209 360 : END SUBROUTINE p86_lda_eval
210 :
211 : ! **************************************************************************************************
212 : !> \brief ...
213 : !> \param rho ...
214 : !> \param rs ...
215 : !> \param grho ...
216 : !> \param e_0 ...
217 : !> \param npoints ...
218 : ! **************************************************************************************************
219 180 : SUBROUTINE p86_u_0(rho, rs, grho, e_0, npoints)
220 :
221 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, rs, grho
222 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
223 : INTEGER, INTENT(in) :: npoints
224 :
225 : INTEGER :: ip
226 : REAL(KIND=dp) :: cr, ep, g, or, phi, r, x
227 :
228 : !$OMP PARALLEL DO PRIVATE(ip,g,r,x,or,cr,phi,ep) DEFAULT(NONE)&
229 180 : !$OMP SHARED(npoints,rho,eps_rho,grho,rs,e_0)
230 : DO ip = 1, npoints
231 :
232 : IF (rho(ip) > eps_rho) THEN
233 : g = grho(ip)
234 : r = rs(ip)
235 : x = r*frs
236 : or = 1.0_dp/rho(ip)
237 : cr = pc1 + (pc2 + a*r + b*r*r)/(1.0_dp + c*r + d*r*r + 1.e4_dp*b*r*r*r)
238 : phi = fpe*pci/cr*g*SQRT(x)*or
239 : ep = EXP(-phi)
240 : e_0(ip) = e_0(ip) + x*or*g*g*cr*ep
241 : END IF
242 :
243 : END DO
244 :
245 180 : END SUBROUTINE p86_u_0
246 :
247 : ! **************************************************************************************************
248 : !> \brief ...
249 : !> \param rho ...
250 : !> \param grho ...
251 : !> \param rs ...
252 : !> \param e_rho ...
253 : !> \param e_ndrho ...
254 : !> \param npoints ...
255 : ! **************************************************************************************************
256 180 : SUBROUTINE p86_u_1(rho, grho, rs, e_rho, e_ndrho, npoints)
257 :
258 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, rs
259 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
260 : INTEGER, INTENT(in) :: npoints
261 :
262 : INTEGER :: ip
263 : REAL(KIND=dp) :: cr, dcr, dphig, dphir, dpv, dq, ep, ff, &
264 : g, or, p, phi, q, r, x
265 :
266 : !$OMP PARALLEL DO PRIVATE(ip,g,r,x,or,p,dpv,q,dq,cr,dcr,dphig,phi,dphir,ep,ff) DEFAULT(NONE)&
267 180 : !$OMP SHARED(npoints,rho,eps_rho,grho,rs,e_rho,e_ndrho)
268 : DO ip = 1, npoints
269 :
270 : IF (rho(ip) > eps_rho) THEN
271 : g = grho(ip)
272 : r = rs(ip)
273 : x = r*frs
274 : or = 1.0_dp/rho(ip)
275 : p = pc2 + a*r + b*r*r
276 : dpv = a + 2.0_dp*b*r
277 : q = 1.0_dp + c*r + d*r*r + 1.e4_dp*b*r*r*r
278 : dq = c + 2.0_dp*d*r + 3.e4_dp*b*r*r
279 : cr = pc1 + p/q
280 : dcr = (dpv*q - p*dq)/(q*q)*(-f13*r*or)
281 : dphig = fpe*pci/cr*SQRT(x)*or
282 : phi = dphig*g
283 : dphir = -phi*(dcr/cr + f76*or)
284 : ep = EXP(-phi)
285 : ff = x*or*g*ep
286 : e_rho(ip) = e_rho(ip) + ff*g*dcr - ff*g*cr*dphir - ff*g*cr*f43*or
287 : e_ndrho(ip) = e_ndrho(ip) + ff*cr*(2.0_dp - g*dphig)
288 : END IF
289 :
290 : END DO
291 :
292 180 : END SUBROUTINE p86_u_1
293 :
294 : ! **************************************************************************************************
295 : !> \brief ...
296 : !> \param rho ...
297 : !> \param grho ...
298 : !> \param rs ...
299 : !> \param e_rho_rho ...
300 : !> \param e_rho_ndrho ...
301 : !> \param e_ndrho_ndrho ...
302 : !> \param npoints ...
303 : ! **************************************************************************************************
304 0 : SUBROUTINE p86_u_2(rho, grho, rs, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
305 : npoints)
306 :
307 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, rs
308 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
309 : INTEGER, INTENT(in) :: npoints
310 :
311 : INTEGER :: ip
312 : REAL(KIND=dp) :: cr, d2cr, d2p, d2phir, d2q, dcr, dphig, &
313 : dphigr, dphir, dpv, dq, ep, g, or, p, &
314 : phi, q, r, x
315 :
316 : !$OMP PARALLEL DO PRIVATE(ip,x,r,cr,phi,ep,g,or,p,q,dpv,dq,dphir,dcr) &
317 : !$OMP PRIVATE(dphig,dphigr,d2phir,d2cr,d2p,d2q) DEFAULT(NONE) &
318 0 : !$OMP SHARED(npoints,rho,eps_rho,grho,rs,e_rho_rho,e_rho_ndrho,e_ndrho_ndrho)
319 : DO ip = 1, npoints
320 :
321 : IF (rho(ip) > eps_rho) THEN
322 : g = grho(ip)
323 : r = rs(ip)
324 : x = r*frs
325 : or = 1.0_dp/rho(ip)
326 : p = pc2 + a*r + b*r*r
327 : dpv = a + 2.0_dp*b*r
328 : d2p = 2.0_dp*b
329 : q = 1.0_dp + c*r + d*r*r + 1.e4_dp*b*r*r*r
330 : dq = c + 2.0_dp*d*r + 3.e4_dp*b*r*r
331 : d2q = 2.0_dp*d + 6.e4_dp*b*r
332 : cr = pc1 + p/q
333 : dcr = (dpv*q - p*dq)/(q*q)*(-f13*r*or)
334 : d2cr = (d2p*q*q - p*q*d2q - 2*dpv*dq*q + 2*p*dq*dq)/(q*q*q)*(f13*r*or)**2 + &
335 : (dpv*q - p*dq)/(q*q)*f13*f43*r*or*or
336 : dphig = fpe*pci/cr*SQRT(x)*or
337 : phi = dphig*g
338 : dphir = -phi*(dcr/cr + f76*or)
339 : d2phir = -dphir*(dcr/cr + f76*or) - &
340 : phi*((d2cr*cr - dcr*dcr)/(cr*cr) - f76*or*or)
341 : dphigr = -dphig*(dcr/cr + f76*or)
342 : ep = EXP(-phi)
343 : e_rho_rho(ip) = e_rho_rho(ip) + x*or*ep*g*g* &
344 : (-f43*or*dcr + d2cr - dcr*dphir + &
345 : f43*or*cr*dphir - dcr*dphir - cr*d2phir + cr*dphir*dphir + &
346 : f43*or*(7.*f13*or*cr - dcr + cr*dphir))
347 : e_rho_ndrho(ip) = e_rho_ndrho(ip) + x*or*ep*g* &
348 : (-2*f43*cr*or + 2*dcr - 2*cr*dphir + f43*or*g*cr*dphig - &
349 : g*dcr*dphig + g*cr*dphir*dphig - g*cr*dphigr)
350 : e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + x*or*ep*cr* &
351 : (2.0_dp - 4.0_dp*g*dphig + g*g*dphig*dphig)
352 : END IF
353 :
354 : END DO
355 :
356 0 : END SUBROUTINE p86_u_2
357 :
358 : ! **************************************************************************************************
359 : !> \brief ...
360 : !> \param rho ...
361 : !> \param grho ...
362 : !> \param rs ...
363 : !> \param e_rho_rho_rho ...
364 : !> \param e_rho_rho_ndrho ...
365 : !> \param e_rho_ndrho_ndrho ...
366 : !> \param e_ndrho_ndrho_ndrho ...
367 : !> \param npoints ...
368 : ! **************************************************************************************************
369 0 : SUBROUTINE p86_u_3(rho, grho, rs, e_rho_rho_rho, &
370 : e_rho_rho_ndrho, e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, &
371 : npoints)
372 :
373 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, rs
374 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
375 : e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
376 : INTEGER, INTENT(in) :: npoints
377 :
378 : INTEGER :: ip
379 : REAL(KIND=dp) :: cr, d2cr, d2p, d2phir, d2phirg, d2pq, d2q, d2z, d3cr, d3phir, d3pq, d3q, &
380 : d3z, dcr, dphig, dphigr, dphir, dpq, dpv, dq, dz, ep, g, or, oz, p, phi, pq, q, r, x
381 :
382 : !$OMP PARALLEL DO PRIVATE(ip,x, r, cr, phi, ep, g, or, p, q, dpv, dq, dphir, dcr, dphig) &
383 : !$OMP PRIVATE(dphigr, d2phir, d3phir, d2cr, d3cr, d2p, d2q, d2phirg, d3q) &
384 : !$OMP PRIVATE(pq, dpq, d2pq, d3pq, oz, dz, d2z, d3z) DEFAULT(NONE) &
385 0 : !$OMP SHARED(npoints,rho,eps_rho,grho,e_rho_rho_rho,e_rho_rho_ndrho,e_rho_ndrho_ndrho,e_ndrho_ndrho_ndrho)
386 : DO ip = 1, npoints
387 :
388 : IF (rho(ip) > eps_rho) THEN
389 : g = grho(ip)
390 : r = rs(ip)
391 : x = r*frs
392 : or = 1.0_dp/rho(ip)
393 : p = pc2 + a*r + b*r*r
394 : dpv = a + 2.0_dp*b*r
395 : d2p = 2.0_dp*b
396 : q = 1.0_dp + c*r + d*r*r + 1.e4_dp*b*r*r*r
397 : dq = c + 2.0_dp*d*r + 3.e4_dp*b*r*r
398 : d2q = 2.0_dp*d + 6.e4*b*r
399 : d3q = 6.e4*b
400 : pq = p/q
401 : dpq = (dpv*q - p*dq)/(q*q)
402 : d2pq = (d2p*q*q - 2*dpv*dq*q + 2*p*dq*dq - p*d2q*q)/(q*q*q)
403 : d3pq = -(3*d2p*dq*q*q - 6*dpv*dq*dq*q + 3*dpv*d2q*q*q + 6*p*dq*dq*dq - 6*p*dq*d2q*q &
404 : + p*d3q*q*q)/(q*q*q*q)
405 : cr = pc1 + pq
406 : dcr = dpq*(-f13*r*or)
407 : d2cr = d2pq*f13*f13*r*r*or*or + dpq*f13*f43*r*or*or
408 : d3cr = d3pq*(-f13*r*or)**3 + 3*d2pq*(-f13*f13*f43*r*r*or*or*or) + &
409 : dpq*(-f13*f43*f13*7*r*or*or*or)
410 : oz = SQRT(x)*or/cr
411 : dz = dcr/cr + f76*or
412 : d2z = d2cr/cr + 2*f76*dcr/cr*or + f76/6.*or*or
413 : d3z = d3cr/cr + 3*f76*d2cr/cr*or + 3*f76/6.*dcr/cr*or*or - 5*f76/36.*or*or*or
414 : dphig = fpe*pci*oz
415 : phi = dphig*g
416 : dphir = -phi*dz
417 : dphigr = -dphig*dz
418 : d2phir = -phi*(d2z - 2*dz*dz)
419 : d3phir = -phi*(d3z - 6*d2z*dz + 6*dz*dz*dz)
420 : d2phirg = -dphigr*dz - &
421 : dphig*((d2cr*cr - dcr*dcr)/(cr*cr) - f76*or*or)
422 : ep = EXP(-phi)
423 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) &
424 : + g*g*x*or*ep*(-280./27.*or*or*or*cr + 3*28./9.*or*or*dcr + &
425 : 3*28./9.*or*or*cr*(-dphir) - 4*or*d2cr - 8*or*dcr*(-dphir) - &
426 : 4*or*cr*(-d2phir + dphir*dphir) + d3cr + 3*d2cr*(-dphir) + &
427 : 3*dcr*(-d2phir + dphir*dphir) + cr*(-d3phir + 3*dphir*d2phir - &
428 : dphir**3))
429 : e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) &
430 : + 2.*x*or*ep*g*(-f43*or*dcr + d2cr - dcr*dphir + &
431 : f43*or*cr*dphir - dcr*dphir - cr*d2phir + cr*dphir*dphir + &
432 : f43*or*(7.*f13*or*cr - dcr + cr*dphir)) - &
433 : dphig*x*or*ep*g*g*(-f43*or*dcr + d2cr - dcr*dphir + &
434 : f43*or*cr*dphir - dcr*dphir - cr*d2phir + cr*dphir*dphir + &
435 : f43*or*(7.*f13*or*cr - dcr + cr*dphir)) + &
436 : x*or*ep*g*g*(-dcr*dphigr + f43*or*cr*dphigr - dcr*dphigr - cr*d2phirg + &
437 : 2.*cr*dphigr*dphir + f43*or*cr*dphigr)
438 : e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) &
439 : + x*or*ep*(-2*f43*cr*or + 2*dcr - 2*cr*dphir + f43*or*g*cr*dphig - &
440 : g*dcr*dphig + g*cr*dphir*dphig - g*cr*dphigr) + &
441 : x*or*ep*g*(-2*cr*dphigr + f43*or*cr*dphig - &
442 : dcr*dphig + cr*dphir*dphig + g*cr*dphigr*dphig - cr*dphigr) - &
443 : x*or*ep*g*dphig*(-2*f43*cr*or + 2*dcr - 2*cr*dphir + f43*or*g*cr*dphig - &
444 : g*dcr*dphig + g*cr*dphir*dphig - g*cr*dphigr)
445 : e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) &
446 : + x*or*ep*cr*dphig*(-6.0_dp + 6.0_dp*g*dphig - g*g*dphig*dphig)
447 : END IF
448 :
449 : END DO
450 :
451 0 : END SUBROUTINE p86_u_3
452 :
453 : END MODULE xc_perdew86
454 :
|
