(git:08825b9)
Loading...
Searching...
No Matches
xc_cs1.F
Go to the documentation of this file.
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 CS1 Functional (Handy s improved LYP functional)
10!> \par History
11!> JGH (26.02.2003) : OpenMP enabled
12!> fawzi (04.2004) : adapted to the new xc interface
13!> \author JGH (16.03.2002)
14! **************************************************************************************************
15MODULE xc_cs1
16
17 USE kinds, ONLY: dp
18 USE xc_derivative_desc, ONLY: deriv_norm_drho,&
19 deriv_norm_drhoa,&
20 deriv_norm_drhob,&
21 deriv_rho,&
22 deriv_rhoa,&
23 deriv_rhob
24 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
25 xc_dset_get_derivative
26 USE xc_derivative_types, ONLY: xc_derivative_get,&
27 xc_derivative_type
28 USE xc_functionals_utilities, ONLY: set_util
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 REAL(KIND=dp), PARAMETER :: f13 = 1.0_dp/3.0_dp, &
39 f23 = 2.0_dp*f13, &
40 f43 = 4.0_dp*f13, &
41 f53 = 5.0_dp*f13
42
43 PUBLIC :: cs1_lda_info, cs1_lda_eval, cs1_lsd_info, cs1_lsd_eval
44
45 REAL(KIND=dp) :: eps_rho
46 LOGICAL :: debug_flag
47 REAL(KIND=dp) :: fsig
48
49 REAL(KIND=dp), PARAMETER :: a = 0.04918_dp, &
50 b = 0.132_dp, &
51 c = 0.2533_dp, &
52 d = 0.349_dp
53
54 REAL(KIND=dp), PARAMETER :: c1 = 0.018897_dp, &
55 c2 = -0.155240_dp, &
56 c3 = 0.159068_dp, &
57 c4 = -0.007953_dp
58
59 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_cs1'
60
61CONTAINS
62
63! **************************************************************************************************
64!> \brief return various information on the functional
65!> \param reference string with the reference of the actual functional
66!> \param shortform string with the shortform of the functional name
67!> \param needs the components needed by this functional are set to
68!> true (does not set the unneeded components to false)
69!> \param max_deriv ...
70! **************************************************************************************************
71 SUBROUTINE cs1_lda_info(reference, shortform, needs, max_deriv)
72 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
73 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
74 INTEGER, INTENT(out), OPTIONAL :: max_deriv
75
76 IF (PRESENT(reference)) THEN
77 reference = "N.C. Handy and A.J. Cohen, J. Chem. Phys., 116, 5411 (2002) {LDA version}"
78 END IF
79 IF (PRESENT(shortform)) THEN
80 shortform = "CS1: Handy improved LYP correlation energy functional {LDA}"
81 END IF
82 IF (PRESENT(needs)) THEN
83 needs%rho = .true.
84 needs%rho_1_3 = .true.
85 needs%norm_drho = .true.
86 END IF
87 IF (PRESENT(max_deriv)) max_deriv = 3
88
89 END SUBROUTINE cs1_lda_info
90
91! **************************************************************************************************
92!> \brief return various information on the functional
93!> \param reference string with the reference of the actual functional
94!> \param shortform string with the shortform of the functional name
95!> \param needs the components needed by this functional are set to
96!> true (does not set the unneeded components to false)
97!> \param max_deriv ...
98! **************************************************************************************************
99 SUBROUTINE cs1_lsd_info(reference, shortform, needs, max_deriv)
100 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
101 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
102 INTEGER, INTENT(out), OPTIONAL :: max_deriv
103
104 IF (PRESENT(reference)) THEN
105 reference = "N.C. Handy and A.J. Cohen, J. Chem. Phys., 116, 5411 (2002)"
106 END IF
107 IF (PRESENT(shortform)) THEN
108 shortform = "CS1: Handy improved LYP correlation energy functional"
109 END IF
110 IF (PRESENT(needs)) THEN
111 needs%rho_spin = .true.
112 needs%rho_spin_1_3 = .true.
113 needs%norm_drho_spin = .true.
114 END IF
115 IF (PRESENT(max_deriv)) max_deriv = 1
116
117 END SUBROUTINE cs1_lsd_info
118
119! **************************************************************************************************
120!> \brief ...
121!> \param cutoff ...
122!> \param debug ...
123! **************************************************************************************************
124 SUBROUTINE cs1_init(cutoff, debug)
125
126 REAL(kind=dp), INTENT(IN) :: cutoff
127 LOGICAL, INTENT(IN), OPTIONAL :: debug
128
129 eps_rho = cutoff
130 CALL set_util(cutoff)
131
132 IF (PRESENT(debug)) THEN
133 debug_flag = debug
134 ELSE
135 debug_flag = .false.
136 END IF
137
138 fsig = 2.0_dp**f13
139
140 END SUBROUTINE cs1_init
141
142! **************************************************************************************************
143!> \brief ...
144!> \param rho_set ...
145!> \param deriv_set ...
146!> \param order ...
147! **************************************************************************************************
148 SUBROUTINE cs1_lda_eval(rho_set, deriv_set, order)
149 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
150 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
151 INTEGER, INTENT(in) :: order
152
153 CHARACTER(len=*), PARAMETER :: routinen = 'cs1_lda_eval'
154
155 INTEGER :: handle, npoints
156 INTEGER, DIMENSION(2, 3) :: bo
157 REAL(kind=dp) :: epsilon_rho
158 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
159 e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
160 e_rho_rho_rho, grho, rho, rho13
161 TYPE(xc_derivative_type), POINTER :: deriv
162
163 CALL timeset(routinen, handle)
164 NULLIFY (rho, rho13, grho, e_0, e_rho, e_ndrho, &
165 e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
166 e_rho_rho_rho, e_rho_rho_ndrho, e_rho_ndrho_ndrho, &
167 e_ndrho_ndrho_ndrho)
168
169 CALL xc_rho_set_get(rho_set, rho_1_3=rho13, rho=rho, &
170 norm_drho=grho, local_bounds=bo, rho_cutoff=epsilon_rho)
171 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
172 CALL cs1_init(epsilon_rho)
173
174 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
175 allocate_deriv=.true.)
176 CALL xc_derivative_get(deriv, deriv_data=e_0)
177
178 CALL cs1_u_0(rho, grho, rho13, e_0, npoints)
179 IF (order >= 1) THEN
180 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
181 allocate_deriv=.true.)
182 CALL xc_derivative_get(deriv, deriv_data=e_rho)
183 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
184 allocate_deriv=.true.)
185 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
186
187 CALL cs1_u_1(rho, grho, rho13, e_rho, e_ndrho, &
188 npoints)
189 END IF
190 IF (order >= 2) THEN
191 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
192 allocate_deriv=.true.)
193 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
194 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
195 allocate_deriv=.true.)
196 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
197 deriv => xc_dset_get_derivative(deriv_set, &
198 [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
199 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
200
201 CALL cs1_u_2(rho, grho, rho13, e_rho_rho, e_rho_ndrho, &
202 e_ndrho_ndrho, npoints)
203 END IF
204 IF (order >= 3) THEN
205 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
206 allocate_deriv=.true.)
207 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
208 deriv => xc_dset_get_derivative(deriv_set, &
209 [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.true.)
210 CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
211 deriv => xc_dset_get_derivative(deriv_set, &
212 [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
213 CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
214 deriv => xc_dset_get_derivative(deriv_set, &
215 [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
216 CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
217
218 CALL cs1_u_3(rho, grho, rho13, e_rho_rho_rho, e_rho_rho_ndrho, &
219 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
220 END IF
221 IF (order > 3) THEN
222 cpabort("derivatives bigger than 3 not implemented")
223 END IF
224
225 CALL timestop(handle)
226 END SUBROUTINE cs1_lda_eval
227
228! **************************************************************************************************
229!> \brief ...
230!> \param rho_set ...
231!> \param deriv_set ...
232!> \param order ...
233! **************************************************************************************************
234 SUBROUTINE cs1_lsd_eval(rho_set, deriv_set, order)
235 TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
236 TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
237 INTEGER, INTENT(in) :: order
238
239 CHARACTER(len=*), PARAMETER :: routinen = 'cs1_lsd_eval'
240
241 INTEGER :: handle, npoints
242 INTEGER, DIMENSION(2, 3) :: bo
243 REAL(kind=dp) :: epsilon_rho
244 REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
245 POINTER :: e_0, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, &
246 norm_drhoa, norm_drhob, rhoa, &
247 rhoa_1_3, rhob, rhob_1_3
248 TYPE(xc_derivative_type), POINTER :: deriv
249
250 CALL timeset(routinen, handle)
251
252 CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
253 rhoa_1_3=rhoa_1_3, rhob_1_3=rhob_1_3, &
254 norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
255 local_bounds=bo, rho_cutoff=epsilon_rho)
256 npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
257 CALL cs1_init(epsilon_rho)
258
259 deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
260 allocate_deriv=.true.)
261 CALL xc_derivative_get(deriv, deriv_data=e_0)
262
263 CALL cs1_ss_0(rhoa, rhob, norm_drhoa, norm_drhob, &
264 rhoa_1_3, rhob_1_3, e_0, npoints)
265 IF (order >= 1) THEN
266 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
267 allocate_deriv=.true.)
268 CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
269 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
270 allocate_deriv=.true.)
271 CALL xc_derivative_get(deriv, deriv_data=e_rhob)
272 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
273 allocate_deriv=.true.)
274 CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
275 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
276 allocate_deriv=.true.)
277 CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
278
279 CALL cs1_ss_1(rhoa, rhob, norm_drhoa, norm_drhob, &
280 rhoa_1_3, rhob_1_3, e_rhoa, e_rhob, e_ndrhoa, e_ndrhob, &
281 npoints)
282 END IF
283 IF (order > 1) THEN
284 cpabort("derivatives bigger than 3 not implemented")
285 END IF
286 CALL timestop(handle)
287 END SUBROUTINE cs1_lsd_eval
288
289! **************************************************************************************************
290!> \brief ...
291!> \param rho ...
292!> \param grho ...
293!> \param r13 ...
294!> \param e_0 ...
295!> \param npoints ...
296! **************************************************************************************************
297 SUBROUTINE cs1_u_0(rho, grho, r13, e_0, npoints)
298
299 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13
300 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
301 INTEGER, INTENT(in) :: npoints
302
303 INTEGER :: ip
304 REAL(kind=dp) :: c2p, c3p, c4p, cp, dpv, f1, f2, f3, f4, &
305 g, oc, ocp, od, odp, r, r3, x
306
307 dpv = fsig*d
308 cp = c*fsig*fsig
309 c2p = c2*fsig**4
310 c3p = -c3*0.25_dp
311 c4p = -c4*0.25_dp
312
313!$OMP PARALLEL DO DEFAULT(NONE) &
314!$OMP SHARED(npoints,rho,eps_rho,e_0,dpv,cp,c2p,c3p,c4p,r13,grho) &
315!$OMP PRIVATE(ip,r,r3,g,x,odp,ocp,od,oc,f1,f2,f3,f4)
316
317 DO ip = 1, npoints
318
319 IF (rho(ip) > eps_rho) THEN
320 r = rho(ip)
321 r3 = r13(ip)
322 g = grho(ip)
323 x = g/(r*r3)
324 odp = 1.0_dp/(r3 + dpv)
325 ocp = 1.0_dp/(r*r*r3*r3 + cp*g*g)
326 od = 1.0_dp/(r3 + d)
327 oc = 1.0_dp/(r*r*r3*r3 + c*g*g)
328 f1 = c1*r*r3*odp
329 f2 = c2p*g**4*r3*r*odp*ocp*ocp
330 f3 = c3p*r*r3*od
331 f4 = c4p*g**4*r3*r*od*oc*oc
332 e_0(ip) = e_0(ip) + f1 + f2 + f3 + f4
333 END IF
334
335 END DO
336
337!$OMP END PARALLEL DO
338
339 END SUBROUTINE cs1_u_0
340
341! **************************************************************************************************
342!> \brief ...
343!> \param rho ...
344!> \param grho ...
345!> \param r13 ...
346!> \param e_rho ...
347!> \param e_ndrho ...
348!> \param npoints ...
349! **************************************************************************************************
350 SUBROUTINE cs1_u_1(rho, grho, r13, e_rho, e_ndrho, npoints)
351
352 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13
353 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
354 INTEGER, INTENT(in) :: npoints
355
356 INTEGER :: ip
357 REAL(kind=dp) :: c2p, c3p, c4p, cp, df1, df3, dgf2, dgf4, &
358 dpv, drf2, drf4, g, oc, ocp, od, odp, &
359 r, r3
360
361 dpv = fsig*d
362 cp = c*fsig*fsig
363 c2p = c2*fsig**4
364 c3p = -c3*0.25_dp
365 c4p = -c4*0.25_dp
366
367!$OMP PARALLEL DO DEFAULT(NONE) &
368!$OMP SHARED(npoints, rho, eps_rho, r13, grho, e_rho) &
369!$OMP SHARED(e_ndrho, dpv, cp, c2p, c3p, c4p) &
370!$OMP PRIVATE(ip, r, r3, g, odp, ocp, od, oc, df1) &
371!$OMP PRIVATE(drf2, dgf2, df3, drf4, dgf4)
372
373 DO ip = 1, npoints
374
375 IF (rho(ip) > eps_rho) THEN
376 r = rho(ip)
377 r3 = r13(ip)
378 g = grho(ip)
379 odp = 1.0_dp/(r3 + dpv)
380 ocp = 1.0_dp/(r*r*r3*r3 + cp*g*g)
381 od = 1.0_dp/(r3 + d)
382 oc = 1.0_dp/(r*r*r3*r3 + c*g*g)
383 df1 = c1*f13*r3*(3*r3 + 4*dpv)*odp*odp
384 drf2 = -f13*c2p*g**4*r3*(13*r**3 - 3*r3*cp*g*g + 12*r*r*r3*r3*dpv - &
385 4*dpv*cp*g*g)*odp**2*ocp**3
386 dgf2 = 4*c2p*g**3*r**4*odp*ocp**3
387 df3 = f13*r3*c3p*(3*r3 + 4*d)*od*od
388 drf4 = -f13*c4p*g**4*r3*(13*r**3 - 3*r3*c*g*g + 12*r*r*r3*r3*d - &
389 4*d*c*g*g)*od**2*oc**3
390 dgf4 = 4*c4p*g**3*r**4*od*oc**3
391 e_rho(ip) = e_rho(ip) + df1 + drf2 + df3 + drf4
392 e_ndrho(ip) = e_ndrho(ip) + dgf2 + dgf4
393 END IF
394
395 END DO
396
397!$OMP END PARALLEL DO
398
399 END SUBROUTINE cs1_u_1
400
401! **************************************************************************************************
402!> \brief ...
403!> \param rho ...
404!> \param grho ...
405!> \param r13 ...
406!> \param e_rho_rho ...
407!> \param e_rho_ndrho ...
408!> \param e_ndrho_ndrho ...
409!> \param npoints ...
410! **************************************************************************************************
411 SUBROUTINE cs1_u_2(rho, grho, r13, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
412 npoints)
413
414 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13
415 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
416 INTEGER, INTENT(in) :: npoints
417
418 INTEGER :: ip
419 REAL(kind=dp) :: c2p, c3p, c4p, cp, d2f1, d2f3, d2gf2, &
420 d2gf4, d2rf2, d2rf4, dpv, drgf2, &
421 drgf4, g, oc, ocp, od, odp, r, r3
422
423 dpv = fsig*d
424 cp = c*fsig*fsig
425 c2p = c2*fsig**4
426 c3p = -c3*0.25_dp
427 c4p = -c4*0.25_dp
428
429!$OMP PARALLEL DO DEFAULT(NONE) &
430!$OMP SHARED(npoints, rho, eps_rho, r13, grho, e_rho_rho) &
431!$OMP SHARED(e_rho_ndrho, e_ndrho_ndrho, dpv, cp, c2p, c3p, c4p) &
432!$OMP PRIVATE(ip,r,r3,g,odp,ocp,od,oc,d2f1,d2rf2,drgf2,d2f3) &
433!$OMP PRIVATE(d2rf4,drgf4,d2gf2,d2gf4)
434
435 DO ip = 1, npoints
436
437 IF (rho(ip) > eps_rho) THEN
438 r = rho(ip)
439 r3 = r13(ip)
440 g = grho(ip)
441 odp = 1.0_dp/(r3 + dpv)
442 ocp = 1.0_dp/(r*r*r3*r3 + cp*g*g)
443 od = 1.0_dp/(r3 + d)
444 oc = 1.0_dp/(r*r*r3*r3 + c*g*g)
445 d2f1 = c1*f23*f13*dpv*r3/r*(r3 + 2*dpv)*odp**3
446 d2rf2 = c2p*f13*f23*g**4*r3/r*(193*dpv*r**5*r3*r3 + 90*dpv*dpv*r**5*r3 &
447 - 88*g*g*cp*r**3*r3 - 100*dpv*dpv*cp*g*g*r*r*r3*r3 &
448 + 2*dpv*dpv*cp*cp*g**4 - 190*g*g*r**3*cp*dpv + g**4*r3*cp*cp*dpv &
449 + 104*r**6)*odp**3*ocp**4
450 drgf2 = c2p*f43*g**3*r*r*r3*(-13*r**3*r3*r3 + 11*cp*r*g*g - 12*dpv*r**3*r3 &
451 + 12*r3*r3*dpv*cp*g*g)*odp*odp*ocp**4
452 d2gf2 = -12*c2p*g*g*r**4*(cp*g*g - r*r*r3*r3)*odp*ocp**4
453 d2f3 = f13*f23*c3p*d*r3/r*(r3 + 2*d)*od**3
454 d2rf4 = c4p*f13*f23*g**4*r3/r*(193*d*r**5*r3*r3 + 90*d*d*r**5*r3 &
455 - 88*g*g*c*r**3*r3 - 100*d*d*c*g*g*r*r*r3*r3 &
456 + 2*d*d*c*c*g**4 - 190*g*g*r**3*c*d + g**4*r3*c*c*d &
457 + 104*r**6)*od**3*oc**4
458 drgf4 = c4p*f43*g**3*r*r*r3*(-13*r**3*r3*r3 + 11*c*r*g*g - 12*d*r**3*r3 &
459 + 12*r3*r3*d*c*g*g)*od*od*oc**4
460 d2gf4 = -12*c4p*g*g*r**4*(c*g*g - r*r*r3*r3)*od*oc**4
461 e_rho_rho(ip) = e_rho_rho(ip) + d2f1 + d2rf2 + d2f3 + d2rf4
462 e_rho_ndrho(ip) = e_rho_ndrho(ip) + drgf2 + drgf4
463 e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + d2gf2 + d2gf4
464 END IF
465
466 END DO
467
468!$OMP END PARALLEL DO
469
470 END SUBROUTINE cs1_u_2
471
472! **************************************************************************************************
473!> \brief ...
474!> \param rho ...
475!> \param grho ...
476!> \param r13 ...
477!> \param e_rho_rho_rho ...
478!> \param e_rho_rho_ndrho ...
479!> \param e_rho_ndrho_ndrho ...
480!> \param e_ndrho_ndrho_ndrho ...
481!> \param npoints ...
482! **************************************************************************************************
483 SUBROUTINE cs1_u_3(rho, grho, r13, e_rho_rho_rho, e_rho_rho_ndrho, &
484 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
485
486 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13
487 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
488 e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
489 INTEGER, INTENT(in) :: npoints
490
491 INTEGER :: ip
492 REAL(kind=dp) :: c2p, c3p, c4p, cp, d2rgf2, d2rgf4, d3f1, d3f3, d3gf2, d3gf4, d3rf2, d3rf4, &
493 dpv, dr2gf2, dr2gf4, g, oc, ocp, od, odp, r, r3, t1, t10, t11, t13, t15, t16, t17, t19, &
494 t2, t22, t23, t26, t29, t3, t32, t33, t37, t4, t44, t45, t50, t51, t52, t58, t6, t61, t7, &
495 t74, t76, t77, t8, t80, t81, t82, t9
496
497 dpv = fsig*d
498 cp = c*fsig*fsig
499 c2p = c2*fsig**4
500 c3p = -c3*0.25_dp
501 c4p = -c4*0.25_dp
502
503!$OMP PARALLEL DO DEFAULT(NONE) &
504!$OMP SHARED(npoints, rho, eps_rho, r13, grho, e_rho_rho_rho) &
505!$OMP SHARED(e_rho_rho_ndrho, e_rho_ndrho_ndrho) &
506!$OMP SHARED(e_ndrho_ndrho_ndrho, dpv, cp, c2p, c3p, c4p) &
507!$OMP PRIVATE(ip,r,r3,g,odp,ocp,od,oc,d3f1,d3rF2,d2rgF2) &
508!$OMP PRIVATE(dr2gF2,d3gF2,d3F3,d3rF4,d2rgF4,dr2gF4,d3gF4) &
509!$OMP PRIVATE(t1, t2, t3, t4, t6, t7, t8, t9, t10, t11, t13) &
510!$OMP PRIVATE(t15, t16, t17, t19, t22, t23, t26, t29, t32) &
511!$OMP PRIVATE(t33, t37, t44, t45, t50, t51, t52, t58, t61) &
512!$OMP PRIVATE(t74, t76, t77, t80, t81, t82)
513
514 DO ip = 1, npoints
515
516 IF (rho(ip) > eps_rho) THEN
517 r = rho(ip)
518 r3 = r13(ip)
519 g = grho(ip)
520 odp = 1.0_dp/(r3 + dpv)
521 ocp = 1.0_dp/(r*r*r3*r3 + cp*g*g)
522 od = 1.0_dp/(r3 + d)
523 oc = 1.0_dp/(r*r*r3*r3 + c*g*g)
524 d3f1 = -c1*f23*f13*f13*dpv*r3/(r*r)*(11*dpv*r3 + 4*dpv*dpv + 4*r/r3)*odp**4
525 t1 = g**2; t2 = t1**2; t3 = r3; t4 = t3**2; t8 = dpv**2; t9 = t8*dpv
526 t10 = cp**2; t11 = t10*cp; t13 = t2*t1; t16 = r**2; t17 = t4*t16
527 t19 = t10*t2; t22 = t16**2; t23 = t22**2; t32 = t22*t16; t37 = t16*r
528 t58 = t22*r; t61 = cp*t1
529 t74 = 4*t9*t11*t13 + 668*t17*t9*t19 + 5524*t4*t23*dpv + 5171*t3*t23*t8 + &
530 1620*t23*t9 - 3728*t3*t32*cp*t1 + 440*t4*t37*t10*t2 + 1500*t2*t3*t37*dpv*t10 &
531 + 4*t13*t4*dpv*t11 + 1737*t37*t8*t19 + 11*t3*t8*t11*t13 - 3860*t3*t58*t9*t61 + &
532 1976*t23*r - 11535*t4*t58*t8*t61 - 11412*t1*t32*cp*dpv
533 t76 = (t3 + dpv)**2; t77 = t76**2; t80 = t17 + t61; t81 = t80**2; t82 = t81**2
534 d3rf2 = -f23*f13*f13*c2p*t2/t4/r*t74/t77/t82/t80
535 t4 = t3*r; t6 = r**2; t7 = t6**2; t8 = t7*t6; t9 = dpv**2; t15 = t1**2
536 t17 = cp**2; t23 = t3**2; t26 = t6*r; t29 = cp*t1; t33 = t17*t15; t44 = t3 + dpv
537 t45 = t44**2; t50 = t23*t6 + t29; t51 = t50**2; t52 = t51**2
538 d2rgf2 = c2p*f23*f43*t1*g*t4*(90*t8*t9 + 193*t3*t8*dpv + 44*t15*t4*t17 - 236*t1 &
539 *t7*cp + 104*t23*t8 - 240*t3*t26*t9*t29 + 54*t23*t9*t33 - 478*t23*t26*dpv*t29 &
540 + 97*r*dpv*t33)/t45/t44/t52/t50
541 dr2gf2 = -4*c2p*g*g*r*r*r3*(-40*r**3*r3*dpv*cp*g*g + 12*r3*r3*dpv*cp*cp*g**4 &
542 + 13*r**6*r3 - 40*r**3*r3*r3*cp*g*g + 11*r*cp*cp*g**4 + 12*r**6*dpv) &
543 *odp*odp*ocp**5
544 d3gf2 = c2p*24*g*r**3*r3*(r**6 - 5*cp*g*g*r**3*r3 + 2*cp*cp*g**4*r3*r3)*odp*ocp**5
545 d3f3 = -f23*f13*f13*c3p*d*r3/(r*r)*(11*d*r3 + 4*d*d + 4*r3*r3)*od**4
546 t1 = g**2; t2 = t1**2; t3 = r3; t4 = t3**2; t8 = d**2; t9 = t8*d
547 t10 = c**2; t11 = t10*c; t13 = t2*t1; t16 = r**2; t17 = t4*t16
548 t19 = t10*t2; t22 = t16**2; t23 = t22**2; t32 = t22*t16; t37 = t16*r
549 t58 = t22*r; t61 = c*t1
550 t74 = 4*t9*t11*t13 + 668*t17*t9*t19 + 5524*t4*t23*d + 5171*t3*t23*t8 + &
551 1620*t23*t9 - 3728*t3*t32*c*t1 + 440*t4*t37*t10*t2 + 1500*t2*t3*t37*d*t10 &
552 + 4*t13*t4*d*t11 + 1737*t37*t8*t19 + 11*t3*t8*t11*t13 - 3860*t3*t58*t9*t61 + &
553 1976*t23*r - 11535*t4*t58*t8*t61 - 11412*t1*t32*c*d
554 t76 = (t3 + d)**2; t77 = t76**2; t80 = t17 + t61; t81 = t80**2; t82 = t81**2
555 d3rf4 = -f23*f13*f13*c4p*t2/t4/r*t74/t77/t82/t80
556 t4 = t3*r; t6 = r**2; t7 = t6**2; t8 = t7*t6; t9 = d**2; t15 = t1**2
557 t17 = c**2; t23 = t3**2; t26 = t6*r; t29 = c*t1; t33 = t17*t15; t44 = t3 + d
558 t45 = t44**2; t50 = t23*t6 + t29; t51 = t50**2; t52 = t51**2
559 d2rgf4 = c4p*f23*f43*t1*g*t4*(90*t8*t9 + 193*t3*t8*d + 44*t15*t4*t17 - 236*t1 &
560 *t7*c + 104*t23*t8 - 240*t3*t26*t9*t29 + 54*t23*t9*t33 - 478*t23*t26*d*t29 &
561 + 97*r*d*t33)/t45/t44/t52/t50
562 dr2gf4 = -4*c4p*g*g*r*r*r3*(-40*r**3*r3*d*c*g*g + 12*r3*r3*d*c*c*g**4 &
563 + 13*r**6*r3 - 40*r**3*r3*r3*c*g*g + 11*r*c*c*g**4 + 12*r**6*d) &
564 *od*od*oc**5
565 d3gf4 = c4p*24*g*r**3*r3*(r**6 - 5*c*g*g*r**3*r3 + 2*c*c*g**4*r3*r3)*od*oc**5
566 e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + d3f1 + d3rf2 + d3f3 + d3rf4
567 e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) + d2rgf2 + d2rgf4
568 e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) + dr2gf2 + dr2gf4
569 e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) + d3gf2 + d3gf4
570 END IF
571
572 END DO
573
574!$OMP END PARALLEL DO
575
576 END SUBROUTINE cs1_u_3
577
578! **************************************************************************************************
579!> \brief ...
580!> \param rhoa ...
581!> \param rhob ...
582!> \param grhoa ...
583!> \param grhob ...
584!> \param r13a ...
585!> \param r13b ...
586!> \param e_0 ...
587!> \param npoints ...
588! **************************************************************************************************
589 SUBROUTINE cs1_ss_0(rhoa, rhob, grhoa, grhob, r13a, r13b, e_0, &
590 npoints)
591
592 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob, grhoa, grhob, r13a, r13b
593 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_0
594 INTEGER, INTENT(in) :: npoints
595
596 INTEGER :: ip
597 REAL(kind=dp) :: f1a, f1b, f2a, f2b, ga, gb, oca, ocb, &
598 oda, odb, r3a, r3b, ra, rb, xa, xb
599
600!FM IF ( .NOT. debug_flag ) CPABORT("Routine not tested")
601
602 cpwarn("not tested!")
603
604!$OMP PARALLEL DO DEFAULT(NONE) &
605!$OMP SHARED(npoints, rhoa, eps_rho, r13a, grhoa, rhob) &
606!$OMP SHARED(r13b, grhob, e_0) &
607!$OMP PRIVATE(ip, f1a, f2a, ra, r3a, ga, xa, oda, oca) &
608!$OMP PRIVATE( f1b, f2b, rb, r3b, gb, xb, odb, ocb)
609
610 DO ip = 1, npoints
611
612 IF (rhoa(ip) < eps_rho) THEN
613 f1a = 0.0_dp
614 f2a = 0.0_dp
615 ELSE
616 ra = rhoa(ip)
617 r3a = r13a(ip)
618 ga = grhoa(ip)
619 xa = ga/(ra*r3a)
620 oda = 1.0_dp/(r3a + d)
621 oca = 1.0_dp/(ra*ra*r3a*r3a + c*ga*ga)
622 f1a = c1*ra*r3a*oda
623 f2a = c2*ga**4*r3a*ra*oda*oca*oca
624 END IF
625 IF (rhob(ip) < eps_rho) THEN
626 f1b = 0.0_dp
627 f2b = 0.0_dp
628 ELSE
629 rb = rhob(ip)
630 r3b = r13b(ip)
631 gb = grhob(ip)
632 xb = gb/(rb*r3b)
633 odb = 1.0_dp/(r3b + d)
634 ocb = 1.0_dp/(rb*rb*r3b*r3b + c*gb*gb)
635 f1b = c1*rb*r3b*odb
636 f2b = c2*gb**4*r3b*rb*odb*ocb*ocb
637 END IF
638
639 e_0(ip) = e_0(ip) + f1a + f1b + f2a + f2b
640
641 END DO
642
643!$OMP END PARALLEL DO
644
645 END SUBROUTINE cs1_ss_0
646
647! **************************************************************************************************
648!> \brief ...
649!> \param rhoa ...
650!> \param rhob ...
651!> \param grhoa ...
652!> \param grhob ...
653!> \param r13a ...
654!> \param r13b ...
655!> \param e_rhoa ...
656!> \param e_rhob ...
657!> \param e_ndrhoa ...
658!> \param e_ndrhob ...
659!> \param npoints ...
660! **************************************************************************************************
661 SUBROUTINE cs1_ss_1(rhoa, rhob, grhoa, grhob, r13a, r13b, e_rhoa, &
662 e_rhob, e_ndrhoa, e_ndrhob, npoints)
663
664 REAL(kind=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob, grhoa, grhob, r13a, r13b
665 REAL(kind=dp), DIMENSION(*), INTENT(INOUT) :: e_rhoa, e_rhob, e_ndrhoa, e_ndrhob
666 INTEGER, INTENT(in) :: npoints
667
668 INTEGER :: ip
669 REAL(kind=dp) :: df1a, df1b, dgf2a, dgf2b, drf2a, drf2b, &
670 ga, gb, oca, ocb, oda, odb, r3a, r3b, &
671 ra, rb
672
673 cpwarn("not tested!")
674
675!$OMP PARALLEL DO DEFAULT(NONE) &
676!$OMP SHARED(npoints, rhoa, eps_rho, r13a, grhoa, rhob) &
677!$OMP SHARED(grhob, e_rhoa, e_ndrhoa, e_rhob, e_ndrhob,r13b) &
678!$OMP PRIVATE(ip, df1a, drf2a, dgf2a, ra, r3a, ga, oda, oca) &
679!$OMP PRIVATE( df1b, drf2b, dgf2b, rb, r3b, gb, odb, ocb)
680
681 DO ip = 1, npoints
682
683 IF (rhoa(ip) < eps_rho) THEN
684 df1a = 0.0_dp
685 drf2a = 0.0_dp
686 dgf2a = 0.0_dp
687 ELSE
688 ra = rhoa(ip)
689 r3a = r13a(ip)
690 ga = grhoa(ip)
691 oda = 1.0_dp/(r3a + d)
692 oca = 1.0_dp/(ra*ra*r3a*r3a + c*ga*ga)
693 df1a = c1*f13*r3a*(3*r3a + 4*d)*oda*oda
694
695 drf2a = -f13*c2*ga**4*r3a*(13*ra**3 - 3*r3a*c*ga*ga + 12*ra*ra*r3a*r3a*d - &
696 4*d*c*ga*ga)*oda**2*oca**3
697
698 dgf2a = 4*c2*ga**3*ra**4*oda*oca**3
699
700 END IF
701 IF (rhob(ip) < eps_rho) THEN
702 df1b = 0.0_dp
703 drf2b = 0.0_dp
704 dgf2b = 0.0_dp
705 ELSE
706 rb = rhob(ip)
707 r3b = r13b(ip)
708 gb = grhob(ip)
709 odb = 1.0_dp/(r3b + d)
710 ocb = 1.0_dp/(rb*rb*r3b*r3b + c*gb*gb)
711 df1b = c1*f13*r3b*(3*r3b + 4*d)*odb*odb
712
713 drf2b = -f13*c2*gb**4*r3b*(13*rb**3 - 3*r3b*c*gb*gb + 12*rb*rb*r3b*r3b*d - &
714 4*d*c*gb*gb)*odb**2*ocb**3
715
716 dgf2b = 4*c2*gb**3*rb**4*odb*ocb**3
717
718 END IF
719
720 e_rhoa(ip) = e_rhoa(ip) + df1a + drf2a
721 e_ndrhoa(ip) = e_ndrhoa(ip) + dgf2a
722 e_rhob(ip) = e_rhob(ip) + df1b + drf2b
723 e_ndrhob(ip) = e_ndrhob(ip) + dgf2b
724
725 END DO
726
727!$OMP END PARALLEL DO
728
729 END SUBROUTINE cs1_ss_1
730
731END MODULE xc_cs1
732
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
xc_cs1
Calculate the CS1 Functional (Handy s improved LYP functional).
Definition xc_cs1.F:15
xc_cs1::cs1_lda_eval
subroutine, public cs1_lda_eval(rho_set, deriv_set, order)
...
Definition xc_cs1.F:149
xc_cs1::cs1_lsd_eval
subroutine, public cs1_lsd_eval(rho_set, deriv_set, order)
...
Definition xc_cs1.F:235
xc_cs1::cs1_lsd_info
subroutine, public cs1_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_cs1.F:100
xc_cs1::cs1_lda_info
subroutine, public cs1_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition xc_cs1.F:72
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_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
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:65
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_functionals_utilities
Utility routines for the functional calculations.
Definition xc_functionals_utilities.F:15
xc_functionals_utilities::set_util
subroutine, public set_util(cutoff)
...
Definition xc_functionals_utilities.F:54
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_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78