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