(git:b279b6b)
xc_cs1.F
Go to the documentation of this file.
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  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, m, 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  m = abs(order)
173  CALL cs1_init(epsilon_rho)
174 
175  IF (order >= 0) THEN
176  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
177  allocate_deriv=.true.)
178  CALL xc_derivative_get(deriv, deriv_data=e_0)
179 
180  CALL cs1_u_0(rho, grho, rho13, e_0, npoints)
181  END IF
182  IF (order >= 1 .OR. order == -1) THEN
183  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
184  allocate_deriv=.true.)
185  CALL xc_derivative_get(deriv, deriv_data=e_rho)
186  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
187  allocate_deriv=.true.)
188  CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
189 
190  CALL cs1_u_1(rho, grho, rho13, e_rho, e_ndrho, &
191  npoints)
192  END IF
193  IF (order >= 2 .OR. order == -2) THEN
194  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho], &
195  allocate_deriv=.true.)
196  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho)
197  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_norm_drho], &
198  allocate_deriv=.true.)
199  CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho)
200  deriv => xc_dset_get_derivative(deriv_set, &
201  [deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
202  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  e_ndrho_ndrho, npoints)
206  END IF
207  IF (order >= 3 .OR. order == -3) THEN
208  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho, deriv_rho, deriv_rho], &
209  allocate_deriv=.true.)
210  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_rho)
211  deriv => xc_dset_get_derivative(deriv_set, &
212  [deriv_rho, deriv_rho, deriv_norm_drho], allocate_deriv=.true.)
213  CALL xc_derivative_get(deriv, deriv_data=e_rho_rho_ndrho)
214  deriv => xc_dset_get_derivative(deriv_set, &
215  [deriv_rho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
216  CALL xc_derivative_get(deriv, deriv_data=e_rho_ndrho_ndrho)
217  deriv => xc_dset_get_derivative(deriv_set, &
218  [deriv_norm_drho, deriv_norm_drho, deriv_norm_drho], allocate_deriv=.true.)
219  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  e_rho_ndrho_ndrho, e_ndrho_ndrho_ndrho, npoints)
223  END IF
224  IF (order > 3 .OR. order < -3) THEN
225  cpabort("derivatives bigger than 3 not implemented")
226  END IF
227 
228  CALL timestop(handle)
229  END SUBROUTINE cs1_lda_eval
230 
231 ! **************************************************************************************************
232 !> \brief ...
233 !> \param rho_set ...
234 !> \param deriv_set ...
235 !> \param order ...
236 ! **************************************************************************************************
237  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  POINTER :: e_0, e_ndrhoa, e_ndrhob, e_rhoa, e_rhob, &
249  norm_drhoa, norm_drhob, rhoa, &
250  rhoa_1_3, rhob, rhob_1_3
251  TYPE(xc_derivative_type), POINTER :: deriv
252 
253  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  local_bounds=bo, rho_cutoff=epsilon_rho)
259  npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
260  CALL cs1_init(epsilon_rho)
261 
262  IF (order >= 0) THEN
263  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
264  allocate_deriv=.true.)
265  CALL xc_derivative_get(deriv, deriv_data=e_0)
266 
267  CALL cs1_ss_0(rhoa, rhob, norm_drhoa, norm_drhob, &
268  rhoa_1_3, rhob_1_3, e_0, npoints)
269  END IF
270  IF (order >= 1 .OR. order == -1) THEN
271  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
272  allocate_deriv=.true.)
273  CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
274  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
275  allocate_deriv=.true.)
276  CALL xc_derivative_get(deriv, deriv_data=e_rhob)
277  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
278  allocate_deriv=.true.)
279  CALL xc_derivative_get(deriv, deriv_data=e_ndrhoa)
280  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
281  allocate_deriv=.true.)
282  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  npoints)
287  END IF
288  IF (order > 1 .OR. order < 1) THEN
289  cpabort("derivatives bigger than 3 not implemented")
290  END IF
291  CALL timestop(handle)
292  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  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  dpv = fsig*d
313  cp = c*fsig*fsig
314  c2p = c2*fsig**4
315  c3p = -c3*0.25_dp
316  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 !$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  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  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  dpv = fsig*d
367  cp = c*fsig*fsig
368  c2p = c2*fsig**4
369  c3p = -c3*0.25_dp
370  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 !$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  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  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  dpv = fsig*d
429  cp = c*fsig*fsig
430  c2p = c2*fsig**4
431  c3p = -c3*0.25_dp
432  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 !$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  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  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  dpv = fsig*d
503  cp = c*fsig*fsig
504  c2p = c2*fsig**4
505  c3p = -c3*0.25_dp
506  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 !$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  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  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  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 !$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  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  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  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 !$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  END SUBROUTINE cs1_ss_1
735 
736 END MODULE xc_cs1
737 
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:238
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:66
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