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 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 : ! **************************************************************************************************
15 : MODULE 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 :
61 : CONTAINS
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 33 : 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 33 : IF (PRESENT(reference)) THEN
77 1 : reference = "N.C. Handy and A.J. Cohen, J. Chem. Phys., 116, 5411 (2002) {LDA version}"
78 : END IF
79 33 : IF (PRESENT(shortform)) THEN
80 1 : shortform = "CS1: Handy improved LYP correlation energy functional {LDA}"
81 : END IF
82 33 : IF (PRESENT(needs)) THEN
83 32 : needs%rho = .TRUE.
84 32 : needs%rho_1_3 = .TRUE.
85 32 : needs%norm_drho = .TRUE.
86 : END IF
87 33 : IF (PRESENT(max_deriv)) max_deriv = 3
88 :
89 33 : 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 0 : 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 0 : IF (PRESENT(reference)) THEN
105 0 : reference = "N.C. Handy and A.J. Cohen, J. Chem. Phys., 116, 5411 (2002)"
106 : END IF
107 0 : IF (PRESENT(shortform)) THEN
108 0 : shortform = "CS1: Handy improved LYP correlation energy functional"
109 : END IF
110 0 : IF (PRESENT(needs)) THEN
111 0 : needs%rho_spin = .TRUE.
112 0 : needs%rho_spin_1_3 = .TRUE.
113 0 : needs%norm_drho_spin = .TRUE.
114 : END IF
115 0 : IF (PRESENT(max_deriv)) max_deriv = 1
116 :
117 0 : END SUBROUTINE cs1_lsd_info
118 :
119 : ! **************************************************************************************************
120 : !> \brief ...
121 : !> \param cutoff ...
122 : !> \param debug ...
123 : ! **************************************************************************************************
124 32 : SUBROUTINE cs1_init(cutoff, debug)
125 :
126 : REAL(KIND=dp), INTENT(IN) :: cutoff
127 : LOGICAL, INTENT(IN), OPTIONAL :: debug
128 :
129 32 : eps_rho = cutoff
130 32 : CALL set_util(cutoff)
131 :
132 32 : IF (PRESENT(debug)) THEN
133 0 : debug_flag = debug
134 : ELSE
135 32 : debug_flag = .FALSE.
136 : END IF
137 :
138 32 : fsig = 2.0_dp**f13
139 :
140 32 : END SUBROUTINE cs1_init
141 :
142 : ! **************************************************************************************************
143 : !> \brief ...
144 : !> \param rho_set ...
145 : !> \param deriv_set ...
146 : !> \param order ...
147 : ! **************************************************************************************************
148 64 : 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, m, npoints
156 : INTEGER, DIMENSION(2, 3) :: bo
157 : REAL(KIND=dp) :: epsilon_rho
158 32 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: e_0, e_ndrho, e_ndrho_ndrho, &
159 32 : e_ndrho_ndrho_ndrho, e_rho, e_rho_ndrho, e_rho_ndrho_ndrho, e_rho_rho, e_rho_rho_ndrho, &
160 32 : e_rho_rho_rho, grho, rho, rho13
161 : TYPE(xc_derivative_type), POINTER :: deriv
162 :
163 32 : CALL timeset(routineN, handle)
164 32 : NULLIFY (rho, rho13, grho, e_0, e_rho, e_ndrho, &
165 32 : e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
166 32 : e_rho_rho_rho, e_rho_rho_ndrho, e_rho_ndrho_ndrho, &
167 32 : e_ndrho_ndrho_ndrho)
168 :
169 : CALL xc_rho_set_get(rho_set, rho_1_3=rho13, rho=rho, &
170 32 : norm_drho=grho, local_bounds=bo, rho_cutoff=epsilon_rho)
171 32 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
172 32 : m = ABS(order)
173 32 : CALL cs1_init(epsilon_rho)
174 :
175 32 : IF (order >= 0) THEN
176 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
177 32 : allocate_deriv=.TRUE.)
178 32 : CALL xc_derivative_get(deriv, deriv_data=e_0)
179 :
180 32 : CALL cs1_u_0(rho, grho, rho13, e_0, npoints)
181 : END IF
182 32 : IF (order >= 1 .OR. order == -1) THEN
183 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
184 32 : allocate_deriv=.TRUE.)
185 32 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
186 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
187 32 : allocate_deriv=.TRUE.)
188 32 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
189 :
190 : CALL cs1_u_1(rho, grho, rho13, e_rho, e_ndrho, &
191 32 : npoints)
192 : END IF
193 32 : IF (order >= 2 .OR. order == -2) THEN
194 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
195 0 : allocate_deriv=.TRUE.)
196 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
197 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
198 0 : allocate_deriv=.TRUE.)
199 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
200 : deriv => xc_dset_get_derivative(deriv_set, &
201 0 : [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
202 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho)
203 :
204 : CALL cs1_u_2(rho, grho, rho13, e_rho_rho, e_rho_ndrho, &
205 0 : e_ndrho_ndrho, npoints)
206 : END IF
207 32 : IF (order >= 3 .OR. order == -3) THEN
208 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
209 0 : allocate_deriv=.TRUE.)
210 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
211 : deriv => xc_dset_get_derivative(deriv_set, &
212 0 : [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.TRUE.)
213 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
214 : deriv => xc_dset_get_derivative(deriv_set, &
215 0 : [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
216 0 : CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
217 : deriv => xc_dset_get_derivative(deriv_set, &
218 0 : [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.TRUE.)
219 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho_ndrho_ndrho)
220 :
221 : CALL cs1_u_3(rho, grho, rho13, e_rho_rho_rho, e_rho_rho_ndrho, &
222 0 : e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
223 : END IF
224 32 : IF (order > 3 .OR. order < -3) THEN
225 0 : CPABORT("derivatives bigger than 3 not implemented")
226 : END IF
227 :
228 32 : CALL timestop(handle)
229 32 : END SUBROUTINE cs1_lda_eval
230 :
231 : ! **************************************************************************************************
232 : !> \brief ...
233 : !> \param rho_set ...
234 : !> \param deriv_set ...
235 : !> \param order ...
236 : ! **************************************************************************************************
237 0 : SUBROUTINE cs1_lsd_eval(rho_set, deriv_set, order)
238 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
239 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
240 : INTEGER, INTENT(in) :: order
241 :
242 : CHARACTER(len=*), PARAMETER :: routineN = 'cs1_lsd_eval'
243 :
244 : INTEGER :: handle, npoints
245 : INTEGER, DIMENSION(2, 3) :: bo
246 : REAL(KIND=dp) :: epsilon_rho
247 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
248 0 : POINTER :: e_0, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, &
249 0 : norm_drhoa, norm_drhob, rhoa, &
250 0 : rhoa_1_3, rhob, rhob_1_3
251 : TYPE(xc_derivative_type), POINTER :: deriv
252 :
253 0 : CALL timeset(routineN, handle)
254 :
255 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
256 : rhoa_1_3=rhoa_1_3, rhob_1_3=rhob_1_3, &
257 : norm_drhoa=norm_drhoa, norm_drhob=norm_drhob, &
258 0 : local_bounds=bo, rho_cutoff=epsilon_rho)
259 0 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
260 0 : CALL cs1_init(epsilon_rho)
261 :
262 0 : IF (order >= 0) THEN
263 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
264 0 : allocate_deriv=.TRUE.)
265 0 : CALL xc_derivative_get(deriv, deriv_data=e_0)
266 :
267 : CALL cs1_ss_0(rhoa, rhob, norm_drhoa, norm_drhob, &
268 0 : rhoa_1_3, rhob_1_3, e_0, npoints)
269 : END IF
270 0 : IF (order >= 1 .OR. order == -1) THEN
271 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
272 0 : allocate_deriv=.TRUE.)
273 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
274 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
275 0 : allocate_deriv=.TRUE.)
276 0 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
277 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
278 0 : allocate_deriv=.TRUE.)
279 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
280 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
281 0 : allocate_deriv=.TRUE.)
282 0 : CALL xc_derivative_get(deriv, deriv_data=e_ndrhob)
283 :
284 : CALL cs1_ss_1(rhoa, rhob, norm_drhoa, norm_drhob, &
285 : rhoa_1_3, rhob_1_3, e_rhoa, e_rhob, e_ndrhoa, e_ndrhob, &
286 0 : npoints)
287 : END IF
288 0 : IF (order > 1 .OR. order < 1) THEN
289 0 : CPABORT("derivatives bigger than 3 not implemented")
290 : END IF
291 0 : CALL timestop(handle)
292 0 : END SUBROUTINE cs1_lsd_eval
293 :
294 : ! **************************************************************************************************
295 : !> \brief ...
296 : !> \param rho ...
297 : !> \param grho ...
298 : !> \param r13 ...
299 : !> \param e_0 ...
300 : !> \param npoints ...
301 : ! **************************************************************************************************
302 32 : SUBROUTINE cs1_u_0(rho, grho, r13, e_0, npoints)
303 :
304 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13
305 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
306 : INTEGER, INTENT(in) :: npoints
307 :
308 : INTEGER :: ip
309 : REAL(KIND=dp) :: c2p, c3p, c4p, cp, dpv, F1, F2, F3, F4, &
310 : g, oc, ocp, od, odp, r, r3, x
311 :
312 32 : dpv = fsig*d
313 32 : cp = c*fsig*fsig
314 32 : c2p = c2*fsig**4
315 32 : c3p = -c3*0.25_dp
316 32 : c4p = -c4*0.25_dp
317 :
318 : !$OMP PARALLEL DO DEFAULT(NONE) &
319 : !$OMP SHARED(npoints,rho,eps_rho,e_0,dpv,cp,c2p,c3p,c4p,r13,grho) &
320 32 : !$OMP PRIVATE(ip,r,r3,g,x,odp,ocp,od,oc,f1,f2,f3,f4)
321 :
322 : DO ip = 1, npoints
323 :
324 : IF (rho(ip) > eps_rho) THEN
325 : r = rho(ip)
326 : r3 = r13(ip)
327 : g = grho(ip)
328 : x = g/(r*r3)
329 : odp = 1.0_dp/(r3 + dpv)
330 : ocp = 1.0_dp/(r*r*r3*r3 + cp*g*g)
331 : od = 1.0_dp/(r3 + d)
332 : oc = 1.0_dp/(r*r*r3*r3 + c*g*g)
333 : F1 = c1*r*r3*odp
334 : F2 = c2p*g**4*r3*r*odp*ocp*ocp
335 : F3 = c3p*r*r3*od
336 : F4 = c4p*g**4*r3*r*od*oc*oc
337 : e_0(ip) = e_0(ip) + F1 + F2 + F3 + F4
338 : END IF
339 :
340 : END DO
341 :
342 : !$OMP END PARALLEL DO
343 :
344 32 : END SUBROUTINE cs1_u_0
345 :
346 : ! **************************************************************************************************
347 : !> \brief ...
348 : !> \param rho ...
349 : !> \param grho ...
350 : !> \param r13 ...
351 : !> \param e_rho ...
352 : !> \param e_ndrho ...
353 : !> \param npoints ...
354 : ! **************************************************************************************************
355 32 : SUBROUTINE cs1_u_1(rho, grho, r13, e_rho, e_ndrho, npoints)
356 :
357 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13
358 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho, e_ndrho
359 : INTEGER, INTENT(in) :: npoints
360 :
361 : INTEGER :: ip
362 : REAL(KIND=dp) :: c2p, c3p, c4p, cp, dF1, dF3, dgF2, dgF4, &
363 : dpv, drF2, drF4, g, oc, ocp, od, odp, &
364 : r, r3
365 :
366 32 : dpv = fsig*d
367 32 : cp = c*fsig*fsig
368 32 : c2p = c2*fsig**4
369 32 : c3p = -c3*0.25_dp
370 32 : c4p = -c4*0.25_dp
371 :
372 : !$OMP PARALLEL DO DEFAULT(NONE) &
373 : !$OMP SHARED(npoints, rho, eps_rho, r13, grho, e_rho) &
374 : !$OMP SHARED(e_ndrho, dpv, cp, c2p, c3p, c4p) &
375 : !$OMP PRIVATE(ip, r, r3, g, odp, ocp, od, oc, df1) &
376 32 : !$OMP PRIVATE(drf2, dgf2, df3, drf4, dgf4)
377 :
378 : DO ip = 1, npoints
379 :
380 : IF (rho(ip) > eps_rho) THEN
381 : r = rho(ip)
382 : r3 = r13(ip)
383 : g = grho(ip)
384 : odp = 1.0_dp/(r3 + dpv)
385 : ocp = 1.0_dp/(r*r*r3*r3 + cp*g*g)
386 : od = 1.0_dp/(r3 + d)
387 : oc = 1.0_dp/(r*r*r3*r3 + c*g*g)
388 : dF1 = c1*f13*r3*(3*r3 + 4*dpv)*odp*odp
389 : drF2 = -f13*c2p*g**4*r3*(13*r**3 - 3*r3*cp*g*g + 12*r*r*r3*r3*dpv - &
390 : 4*dpv*cp*g*g)*odp**2*ocp**3
391 : dgF2 = 4*c2p*g**3*r**4*odp*ocp**3
392 : dF3 = f13*r3*c3p*(3*r3 + 4*d)*od*od
393 : drF4 = -f13*c4p*g**4*r3*(13*r**3 - 3*r3*c*g*g + 12*r*r*r3*r3*d - &
394 : 4*d*c*g*g)*od**2*oc**3
395 : dgF4 = 4*c4p*g**3*r**4*od*oc**3
396 : e_rho(ip) = e_rho(ip) + dF1 + drF2 + dF3 + drF4
397 : e_ndrho(ip) = e_ndrho(ip) + dgF2 + dgF4
398 : END IF
399 :
400 : END DO
401 :
402 : !$OMP END PARALLEL DO
403 :
404 32 : END SUBROUTINE cs1_u_1
405 :
406 : ! **************************************************************************************************
407 : !> \brief ...
408 : !> \param rho ...
409 : !> \param grho ...
410 : !> \param r13 ...
411 : !> \param e_rho_rho ...
412 : !> \param e_rho_ndrho ...
413 : !> \param e_ndrho_ndrho ...
414 : !> \param npoints ...
415 : ! **************************************************************************************************
416 0 : SUBROUTINE cs1_u_2(rho, grho, r13, e_rho_rho, e_rho_ndrho, e_ndrho_ndrho, &
417 : npoints)
418 :
419 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13
420 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho, e_rho_ndrho, e_ndrho_ndrho
421 : INTEGER, INTENT(in) :: npoints
422 :
423 : INTEGER :: ip
424 : REAL(KIND=dp) :: c2p, c3p, c4p, cp, d2F1, d2F3, d2gF2, &
425 : d2gF4, d2rF2, d2rF4, dpv, drgF2, &
426 : drgF4, g, oc, ocp, od, odp, r, r3
427 :
428 0 : dpv = fsig*d
429 0 : cp = c*fsig*fsig
430 0 : c2p = c2*fsig**4
431 0 : c3p = -c3*0.25_dp
432 0 : c4p = -c4*0.25_dp
433 :
434 : !$OMP PARALLEL DO DEFAULT(NONE) &
435 : !$OMP SHARED(npoints, rho, eps_rho, r13, grho, e_rho_rho) &
436 : !$OMP SHARED(e_rho_ndrho, e_ndrho_ndrho, dpv, cp, c2p, c3p, c4p) &
437 : !$OMP PRIVATE(ip,r,r3,g,odp,ocp,od,oc,d2f1,d2rf2,drgf2,d2f3) &
438 0 : !$OMP PRIVATE(d2rf4,drgf4,d2gf2,d2gf4)
439 :
440 : DO ip = 1, npoints
441 :
442 : IF (rho(ip) > eps_rho) THEN
443 : r = rho(ip)
444 : r3 = r13(ip)
445 : g = grho(ip)
446 : odp = 1.0_dp/(r3 + dpv)
447 : ocp = 1.0_dp/(r*r*r3*r3 + cp*g*g)
448 : od = 1.0_dp/(r3 + d)
449 : oc = 1.0_dp/(r*r*r3*r3 + c*g*g)
450 : d2F1 = c1*f23*f13*dpv*r3/r*(r3 + 2*dpv)*odp**3
451 : d2rF2 = c2p*f13*f23*g**4*r3/r*(193*dpv*r**5*r3*r3 + 90*dpv*dpv*r**5*r3 &
452 : - 88*g*g*cp*r**3*r3 - 100*dpv*dpv*cp*g*g*r*r*r3*r3 &
453 : + 2*dpv*dpv*cp*cp*g**4 - 190*g*g*r**3*cp*dpv + g**4*r3*cp*cp*dpv &
454 : + 104*r**6)*odp**3*ocp**4
455 : drgF2 = c2p*f43*g**3*r*r*r3*(-13*r**3*r3*r3 + 11*cp*r*g*g - 12*dpv*r**3*r3 &
456 : + 12*r3*r3*dpv*cp*g*g)*odp*odp*ocp**4
457 : d2gF2 = -12*c2p*g*g*r**4*(cp*g*g - r*r*r3*r3)*odp*ocp**4
458 : d2F3 = f13*f23*c3p*d*r3/r*(r3 + 2*d)*od**3
459 : d2rF4 = c4p*f13*f23*g**4*r3/r*(193*d*r**5*r3*r3 + 90*d*d*r**5*r3 &
460 : - 88*g*g*c*r**3*r3 - 100*d*d*c*g*g*r*r*r3*r3 &
461 : + 2*d*d*c*c*g**4 - 190*g*g*r**3*c*d + g**4*r3*c*c*d &
462 : + 104*r**6)*od**3*oc**4
463 : drgF4 = c4p*f43*g**3*r*r*r3*(-13*r**3*r3*r3 + 11*c*r*g*g - 12*d*r**3*r3 &
464 : + 12*r3*r3*d*c*g*g)*od*od*oc**4
465 : d2gF4 = -12*c4p*g*g*r**4*(c*g*g - r*r*r3*r3)*od*oc**4
466 : e_rho_rho(ip) = e_rho_rho(ip) + d2F1 + d2rF2 + d2F3 + d2rF4
467 : e_rho_ndrho(ip) = e_rho_ndrho(ip) + drgF2 + drgF4
468 : e_ndrho_ndrho(ip) = e_ndrho_ndrho(ip) + d2gF2 + d2gF4
469 : END IF
470 :
471 : END DO
472 :
473 : !$OMP END PARALLEL DO
474 :
475 0 : END SUBROUTINE cs1_u_2
476 :
477 : ! **************************************************************************************************
478 : !> \brief ...
479 : !> \param rho ...
480 : !> \param grho ...
481 : !> \param r13 ...
482 : !> \param e_rho_rho_rho ...
483 : !> \param e_rho_rho_ndrho ...
484 : !> \param e_rho_ndrho_ndrho ...
485 : !> \param e_ndrho_ndrho_ndrho ...
486 : !> \param npoints ...
487 : ! **************************************************************************************************
488 0 : SUBROUTINE cs1_u_3(rho, grho, r13, e_rho_rho_rho, e_rho_rho_ndrho, &
489 : e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
490 :
491 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, grho, r13
492 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rho_rho_rho, e_rho_rho_ndrho, &
493 : e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho
494 : INTEGER, INTENT(in) :: npoints
495 :
496 : INTEGER :: ip
497 : REAL(KIND=dp) :: c2p, c3p, c4p, cp, d2rgF2, d2rgF4, d3F1, d3F3, d3gF2, d3gF4, d3rF2, d3rF4, &
498 : dpv, dr2gF2, dr2gF4, g, oc, ocp, od, odp, r, r3, t1, t10, t11, t13, t15, t16, t17, t19, &
499 : t2, t22, t23, t26, t29, t3, t32, t33, t37, t4, t44, t45, t50, t51, t52, t58, t6, t61, t7, &
500 : t74, t76, t77, t8, t80, t81, t82, t9
501 :
502 0 : dpv = fsig*d
503 0 : cp = c*fsig*fsig
504 0 : c2p = c2*fsig**4
505 0 : c3p = -c3*0.25_dp
506 0 : c4p = -c4*0.25_dp
507 :
508 : !$OMP PARALLEL DO DEFAULT(NONE) &
509 : !$OMP SHARED(npoints, rho, eps_rho, r13, grho, e_rho_rho_rho) &
510 : !$OMP SHARED(e_rho_rho_ndrho, e_rho_ndrho_ndrho) &
511 : !$OMP SHARED(e_ndrho_ndrho_ndrho, dpv, cp, c2p, c3p, c4p) &
512 : !$OMP PRIVATE(ip,r,r3,g,odp,ocp,od,oc,d3f1,d3rF2,d2rgF2) &
513 : !$OMP PRIVATE(dr2gF2,d3gF2,d3F3,d3rF4,d2rgF4,dr2gF4,d3gF4) &
514 : !$OMP PRIVATE(t1, t2, t3, t4, t6, t7, t8, t9, t10, t11, t13) &
515 : !$OMP PRIVATE(t15, t16, t17, t19, t22, t23, t26, t29, t32) &
516 : !$OMP PRIVATE(t33, t37, t44, t45, t50, t51, t52, t58, t61) &
517 0 : !$OMP PRIVATE(t74, t76, t77, t80, t81, t82)
518 :
519 : DO ip = 1, npoints
520 :
521 : IF (rho(ip) > eps_rho) THEN
522 : r = rho(ip)
523 : r3 = r13(ip)
524 : g = grho(ip)
525 : odp = 1.0_dp/(r3 + dpv)
526 : ocp = 1.0_dp/(r*r*r3*r3 + cp*g*g)
527 : od = 1.0_dp/(r3 + d)
528 : oc = 1.0_dp/(r*r*r3*r3 + c*g*g)
529 : d3F1 = -c1*f23*f13*f13*dpv*r3/(r*r)*(11*dpv*r3 + 4*dpv*dpv + 4*r/r3)*odp**4
530 : t1 = g**2; t2 = t1**2; t3 = r3; t4 = t3**2; t8 = dpv**2; t9 = t8*dpv
531 : t10 = cp**2; t11 = t10*cp; t13 = t2*t1; t16 = r**2; t17 = t4*t16
532 : t19 = t10*t2; t22 = t16**2; t23 = t22**2; t32 = t22*t16; t37 = t16*r
533 : t58 = t22*r; t61 = cp*t1
534 : t74 = 4*t9*t11*t13 + 668*t17*t9*t19 + 5524*t4*t23*dpv + 5171*t3*t23*t8 + &
535 : 1620*t23*t9 - 3728*t3*t32*cp*t1 + 440*t4*t37*t10*t2 + 1500*t2*t3*t37*dpv*t10 &
536 : + 4*t13*t4*dpv*t11 + 1737*t37*t8*t19 + 11*t3*t8*t11*t13 - 3860*t3*t58*t9*t61 + &
537 : 1976*t23*r - 11535*t4*t58*t8*t61 - 11412*t1*t32*cp*dpv
538 : t76 = (t3 + dpv)**2; t77 = t76**2; t80 = t17 + t61; t81 = t80**2; t82 = t81**2
539 : d3rF2 = -f23*f13*f13*c2p*t2/t4/r*t74/t77/t82/t80
540 : t4 = t3*r; t6 = r**2; t7 = t6**2; t8 = t7*t6; t9 = dpv**2; t15 = t1**2
541 : t17 = cp**2; t23 = t3**2; t26 = t6*r; t29 = cp*t1; t33 = t17*t15; t44 = t3 + dpv
542 : t45 = t44**2; t50 = t23*t6 + t29; t51 = t50**2; t52 = t51**2
543 : d2rgF2 = c2p*f23*f43*t1*g*t4*(90*t8*t9 + 193*t3*t8*dpv + 44*t15*t4*t17 - 236*t1 &
544 : *t7*cp + 104*t23*t8 - 240*t3*t26*t9*t29 + 54*t23*t9*t33 - 478*t23*t26*dpv*t29 &
545 : + 97*r*dpv*t33)/t45/t44/t52/t50
546 : dr2gF2 = -4*c2p*g*g*r*r*r3*(-40*r**3*r3*dpv*cp*g*g + 12*r3*r3*dpv*cp*cp*g**4 &
547 : + 13*r**6*r3 - 40*r**3*r3*r3*cp*g*g + 11*r*cp*cp*g**4 + 12*r**6*dpv) &
548 : *odp*odp*ocp**5
549 : 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
550 : d3F3 = -f23*f13*f13*c3p*d*r3/(r*r)*(11*d*r3 + 4*d*d + 4*r3*r3)*od**4
551 : t1 = g**2; t2 = t1**2; t3 = r3; t4 = t3**2; t8 = d**2; t9 = t8*d
552 : t10 = c**2; t11 = t10*c; t13 = t2*t1; t16 = r**2; t17 = t4*t16
553 : t19 = t10*t2; t22 = t16**2; t23 = t22**2; t32 = t22*t16; t37 = t16*r
554 : t58 = t22*r; t61 = c*t1
555 : t74 = 4*t9*t11*t13 + 668*t17*t9*t19 + 5524*t4*t23*d + 5171*t3*t23*t8 + &
556 : 1620*t23*t9 - 3728*t3*t32*c*t1 + 440*t4*t37*t10*t2 + 1500*t2*t3*t37*d*t10 &
557 : + 4*t13*t4*d*t11 + 1737*t37*t8*t19 + 11*t3*t8*t11*t13 - 3860*t3*t58*t9*t61 + &
558 : 1976*t23*r - 11535*t4*t58*t8*t61 - 11412*t1*t32*c*d
559 : t76 = (t3 + d)**2; t77 = t76**2; t80 = t17 + t61; t81 = t80**2; t82 = t81**2
560 : d3rF4 = -f23*f13*f13*c4p*t2/t4/r*t74/t77/t82/t80
561 : t4 = t3*r; t6 = r**2; t7 = t6**2; t8 = t7*t6; t9 = d**2; t15 = t1**2
562 : t17 = c**2; t23 = t3**2; t26 = t6*r; t29 = c*t1; t33 = t17*t15; t44 = t3 + d
563 : t45 = t44**2; t50 = t23*t6 + t29; t51 = t50**2; t52 = t51**2
564 : d2rgF4 = c4p*f23*f43*t1*g*t4*(90*t8*t9 + 193*t3*t8*d + 44*t15*t4*t17 - 236*t1 &
565 : *t7*c + 104*t23*t8 - 240*t3*t26*t9*t29 + 54*t23*t9*t33 - 478*t23*t26*d*t29 &
566 : + 97*r*d*t33)/t45/t44/t52/t50
567 : dr2gF4 = -4*c4p*g*g*r*r*r3*(-40*r**3*r3*d*c*g*g + 12*r3*r3*d*c*c*g**4 &
568 : + 13*r**6*r3 - 40*r**3*r3*r3*c*g*g + 11*r*c*c*g**4 + 12*r**6*d) &
569 : *od*od*oc**5
570 : 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
571 : e_rho_rho_rho(ip) = e_rho_rho_rho(ip) + d3F1 + d3rF2 + d3F3 + d3rF4
572 : e_rho_rho_ndrho(ip) = e_rho_rho_ndrho(ip) + d2rgF2 + d2rgF4
573 : e_rho_ndrho_ndrho(ip) = e_rho_ndrho_ndrho(ip) + dr2gF2 + dr2gF4
574 : e_ndrho_ndrho_ndrho(ip) = e_ndrho_ndrho_ndrho(ip) + d3gF2 + d3gF4
575 : END IF
576 :
577 : END DO
578 :
579 : !$OMP END PARALLEL DO
580 :
581 0 : END SUBROUTINE cs1_u_3
582 :
583 : ! **************************************************************************************************
584 : !> \brief ...
585 : !> \param rhoa ...
586 : !> \param rhob ...
587 : !> \param grhoa ...
588 : !> \param grhob ...
589 : !> \param r13a ...
590 : !> \param r13b ...
591 : !> \param e_0 ...
592 : !> \param npoints ...
593 : ! **************************************************************************************************
594 0 : SUBROUTINE cs1_ss_0(rhoa, rhob, grhoa, grhob, r13a, r13b, e_0, &
595 : npoints)
596 :
597 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob, grhoa, grhob, r13a, r13b
598 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0
599 : INTEGER, INTENT(in) :: npoints
600 :
601 : INTEGER :: ip
602 : REAL(KIND=dp) :: F1a, F1b, F2a, F2b, ga, gb, oca, ocb, &
603 : oda, odb, r3a, r3b, ra, rb, xa, xb
604 :
605 : !FM IF ( .NOT. debug_flag ) CPABORT("Routine not tested")
606 :
607 0 : CPWARN("not tested!")
608 :
609 : !$OMP PARALLEL DO DEFAULT(NONE) &
610 : !$OMP SHARED(npoints, rhoa, eps_rho, r13a, grhoa, rhob) &
611 : !$OMP SHARED(r13b, grhob, e_0) &
612 : !$OMP PRIVATE(ip, f1a, f2a, ra, r3a, ga, xa, oda, oca) &
613 0 : !$OMP PRIVATE( f1b, f2b, rb, r3b, gb, xb, odb, ocb)
614 :
615 : DO ip = 1, npoints
616 :
617 : IF (rhoa(ip) < eps_rho) THEN
618 : F1a = 0.0_dp
619 : F2a = 0.0_dp
620 : ELSE
621 : ra = rhoa(ip)
622 : r3a = r13a(ip)
623 : ga = grhoa(ip)
624 : xa = ga/(ra*r3a)
625 : oda = 1.0_dp/(r3a + d)
626 : oca = 1.0_dp/(ra*ra*r3a*r3a + c*ga*ga)
627 : F1a = c1*ra*r3a*oda
628 : F2a = c2*ga**4*r3a*ra*oda*oca*oca
629 : END IF
630 : IF (rhob(ip) < eps_rho) THEN
631 : F1b = 0.0_dp
632 : F2b = 0.0_dp
633 : ELSE
634 : rb = rhob(ip)
635 : r3b = r13b(ip)
636 : gb = grhob(ip)
637 : xb = gb/(rb*r3b)
638 : odb = 1.0_dp/(r3b + d)
639 : ocb = 1.0_dp/(rb*rb*r3b*r3b + c*gb*gb)
640 : F1b = c1*rb*r3b*odb
641 : F2b = c2*gb**4*r3b*rb*odb*ocb*ocb
642 : END IF
643 :
644 : e_0(ip) = e_0(ip) + F1a + F1b + F2a + F2b
645 :
646 : END DO
647 :
648 : !$OMP END PARALLEL DO
649 :
650 0 : END SUBROUTINE cs1_ss_0
651 :
652 : ! **************************************************************************************************
653 : !> \brief ...
654 : !> \param rhoa ...
655 : !> \param rhob ...
656 : !> \param grhoa ...
657 : !> \param grhob ...
658 : !> \param r13a ...
659 : !> \param r13b ...
660 : !> \param e_rhoa ...
661 : !> \param e_rhob ...
662 : !> \param e_ndrhoa ...
663 : !> \param e_ndrhob ...
664 : !> \param npoints ...
665 : ! **************************************************************************************************
666 0 : SUBROUTINE cs1_ss_1(rhoa, rhob, grhoa, grhob, r13a, r13b, e_rhoa, &
667 : e_rhob, e_ndrhoa, e_ndrhob, npoints)
668 :
669 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rhoa, rhob, grhoa, grhob, r13a, r13b
670 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_rhoa, e_rhob, e_ndrhoa, e_ndrhob
671 : INTEGER, INTENT(in) :: npoints
672 :
673 : INTEGER :: ip
674 : REAL(KIND=dp) :: dF1a, dF1b, dgF2a, dgF2b, drF2a, drF2b, &
675 : ga, gb, oca, ocb, oda, odb, r3a, r3b, &
676 : ra, rb
677 :
678 0 : CPWARN("not tested!")
679 :
680 : !$OMP PARALLEL DO DEFAULT(NONE) &
681 : !$OMP SHARED(npoints, rhoa, eps_rho, r13a, grhoa, rhob) &
682 : !$OMP SHARED(grhob, e_rhoa, e_ndrhoa, e_rhob, e_ndrhob,r13b) &
683 : !$OMP PRIVATE(ip, df1a, drf2a, dgf2a, ra, r3a, ga, oda, oca) &
684 0 : !$OMP PRIVATE( df1b, drf2b, dgf2b, rb, r3b, gb, odb, ocb)
685 :
686 : DO ip = 1, npoints
687 :
688 : IF (rhoa(ip) < eps_rho) THEN
689 : dF1a = 0.0_dp
690 : drF2a = 0.0_dp
691 : dgF2a = 0.0_dp
692 : ELSE
693 : ra = rhoa(ip)
694 : r3a = r13a(ip)
695 : ga = grhoa(ip)
696 : oda = 1.0_dp/(r3a + d)
697 : oca = 1.0_dp/(ra*ra*r3a*r3a + c*ga*ga)
698 : dF1a = c1*f13*r3a*(3*r3a + 4*d)*oda*oda
699 :
700 : drF2a = -f13*c2*ga**4*r3a*(13*ra**3 - 3*r3a*c*ga*ga + 12*ra*ra*r3a*r3a*d - &
701 : 4*d*c*ga*ga)*oda**2*oca**3
702 :
703 : dgF2a = 4*c2*ga**3*ra**4*oda*oca**3
704 :
705 : END IF
706 : IF (rhob(ip) < eps_rho) THEN
707 : dF1b = 0.0_dp
708 : drF2b = 0.0_dp
709 : dgF2b = 0.0_dp
710 : ELSE
711 : rb = rhob(ip)
712 : r3b = r13b(ip)
713 : gb = grhob(ip)
714 : odb = 1.0_dp/(r3b + d)
715 : ocb = 1.0_dp/(rb*rb*r3b*r3b + c*gb*gb)
716 : dF1b = c1*f13*r3b*(3*r3b + 4*d)*odb*odb
717 :
718 : drF2b = -f13*c2*gb**4*r3b*(13*rb**3 - 3*r3b*c*gb*gb + 12*rb*rb*r3b*r3b*d - &
719 : 4*d*c*gb*gb)*odb**2*ocb**3
720 :
721 : dgF2b = 4*c2*gb**3*rb**4*odb*ocb**3
722 :
723 : END IF
724 :
725 : e_rhoa(ip) = e_rhoa(ip) + dF1a + drF2a
726 : e_ndrhoa(ip) = e_ndrhoa(ip) + dgF2a
727 : e_rhob(ip) = e_rhob(ip) + dF1b + drF2b
728 : e_ndrhob(ip) = e_ndrhob(ip) + dgF2b
729 :
730 : END DO
731 :
732 : !$OMP END PARALLEL DO
733 :
734 0 : END SUBROUTINE cs1_ss_1
735 :
736 : END MODULE xc_cs1
737 :
|
