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 Hamprecht, Cohen, Tozer, and Handy (HCTH) exchange
10 : !> functional
11 : !> \author fawzi
12 : ! **************************************************************************************************
13 : MODULE xc_hcth
14 :
15 : USE cp_log_handling, ONLY: cp_to_string
16 : USE kinds, ONLY: dp
17 : USE mathconstants, ONLY: pi
18 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
19 : deriv_rho
20 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
21 : xc_dset_get_derivative
22 : USE xc_derivative_types, ONLY: xc_derivative_get,&
23 : xc_derivative_type
24 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
25 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
26 : xc_rho_set_type
27 : #include "../base/base_uses.f90"
28 :
29 : IMPLICIT NONE
30 : PRIVATE
31 :
32 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_hcth'
34 :
35 : PUBLIC :: hcth_lda_info, hcth_lda_eval
36 : CONTAINS
37 :
38 : ! **************************************************************************************************
39 : !> \brief return various information on the functional
40 : !> \param iparset ...
41 : !> \param reference string with the reference of the actual functional
42 : !> \param shortform string with the shortform of the functional name
43 : !> \param needs the components needed by this functional are set to
44 : !> true (does not set the unneeded components to false)
45 : !> \param max_deriv ...
46 : !> \author fawzi
47 : ! **************************************************************************************************
48 589 : SUBROUTINE hcth_lda_info(iparset, reference, shortform, needs, max_deriv)
49 : INTEGER, INTENT(in) :: iparset
50 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
51 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
52 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
53 :
54 589 : SELECT CASE (iparset)
55 : CASE (93)
56 0 : IF (PRESENT(reference)) THEN
57 : reference = "F. A. Hamprecht, A. J. Cohen, D. J. Tozer, and N. C. Handy, J. Chem. Phys. 109, 6264 (1998);"// &
58 0 : " HCTH/93 xc functional {LDA version}"
59 : END IF
60 0 : IF (PRESENT(shortform)) THEN
61 0 : shortform = "HCTH/93 xc energy functional (LDA)"
62 : END IF
63 : CASE (120)
64 589 : IF (PRESENT(reference)) THEN
65 : reference = "A. D. Boese, N. L. Doltsinis, N. C. Handy, and M. Sprik, J. Chem. Phys. 112, 1670 (2000);"// &
66 2 : " HCTH/120 xc functional {LDA version}"
67 : END IF
68 589 : IF (PRESENT(shortform)) THEN
69 2 : shortform = "HCTH/120 xc energy functional (LDA)"
70 : END IF
71 : CASE (147)
72 0 : IF (PRESENT(reference)) THEN
73 : reference = "A. D. Boese, N. L. Doltsinis, N. C. Handy, and M. Sprik, J. Chem. Phys. 112, 1670 (2000);"// &
74 0 : " HCTH/147 xc functional {LDA Version}"
75 : END IF
76 0 : IF (PRESENT(shortform)) THEN
77 0 : shortform = "HCTH/147 xc energy functional (LDA)"
78 : END IF
79 : CASE (407)
80 0 : IF (PRESENT(reference)) THEN
81 : reference = "A. D. Boese and N. C. Handy, J. Chem. Phys. 114, 5497 (2001); "// &
82 0 : "HCTH/407 xc functional {LDA version}"
83 : END IF
84 0 : IF (PRESENT(shortform)) THEN
85 0 : shortform = "HCTH/407 xc energy functional (LDA)"
86 : END IF
87 : CASE (408)
88 0 : IF (PRESENT(reference)) THEN
89 : reference = "P. Verma and D. G. Truhlar, J. Phys. Chem. Lett. 8, 380 (2016); "// &
90 0 : "HLE16 xc functional {LDA version}"
91 : END IF
92 0 : IF (PRESENT(shortform)) THEN
93 0 : shortform = "HLE16 xc energy functional (LDA)"
94 : END IF
95 : CASE default
96 0 : CPABORT("Invalid HCTH parameter set requested ("//cp_to_string(iparset)//")")
97 : END SELECT
98 589 : IF (PRESENT(needs)) THEN
99 587 : needs%rho = .TRUE.
100 587 : needs%norm_drho = .TRUE.
101 : END IF
102 589 : IF (PRESENT(max_deriv)) max_deriv = 1
103 :
104 589 : END SUBROUTINE hcth_lda_info
105 :
106 : ! **************************************************************************************************
107 : !> \brief evaluates the hcth functional for lda
108 : !> \param iparset the parameter set that should be used (93,120,147,407)
109 : !> \param rho_set the density where you want to evaluate the functional
110 : !> \param deriv_set place where to store the functional derivatives (they are
111 : !> added to the derivatives)
112 : !> \param grad_deriv degree of the derivative that should be evaluated,
113 : !> if positive all the derivatives up to the given degree are evaluated,
114 : !> if negative only the given degree is calculated
115 : !> \author fawzi
116 : ! **************************************************************************************************
117 565 : SUBROUTINE hcth_lda_eval(iparset, rho_set, deriv_set, grad_deriv)
118 : INTEGER, INTENT(in) :: iparset
119 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
120 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
121 : INTEGER, INTENT(in) :: grad_deriv
122 :
123 : INTEGER :: npoints
124 : INTEGER, DIMENSION(2, 3) :: bo
125 : REAL(kind=dp) :: epsilon_rho
126 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
127 565 : POINTER :: e_0, e_ndrho, e_rho, norm_drho, rho
128 : TYPE(xc_derivative_type), POINTER :: deriv
129 :
130 565 : NULLIFY (e_0, e_ndrho, e_rho, norm_drho, rho)
131 :
132 : CALL xc_rho_set_get(rho_set, rho=rho, &
133 565 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
134 565 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
135 :
136 565 : IF (grad_deriv >= 0) THEN
137 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
138 565 : allocate_deriv=.TRUE.)
139 565 : CALL xc_derivative_get(deriv, deriv_data=e_0)
140 : END IF
141 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
142 565 : allocate_deriv=.TRUE.)
143 565 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
144 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
145 565 : allocate_deriv=.TRUE.)
146 565 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
147 565 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
148 0 : CPABORT("derivatives bigger than 1 not implemented")
149 : END IF
150 :
151 : CALL hcth_lda_calc(iparset=iparset, rho=rho, norm_drho=norm_drho, &
152 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
153 565 : npoints=npoints, epsilon_rho=epsilon_rho)
154 565 : END SUBROUTINE hcth_lda_eval
155 :
156 : ! **************************************************************************************************
157 : !> \brief Calculate the gradient-corrected xc energy and potential
158 : !> of Hamprecht, Cohen, Tozer, and Handy (HCTH) for a closed shell
159 : !> density.
160 : !> \param iparset the parameter set that should be used (93,120,147,407)
161 : !> \param rho the density
162 : !> \param norm_drho the norm of the gradient of the density
163 : !> \param e_0 the value of the functional in that point
164 : !> \param e_rho the derivative of the functional wrt. rho
165 : !> \param e_ndrho the derivative of the functional wrt. norm_drho
166 : !> \param epsilon_rho the cutoff on rho
167 : !> \param npoints ...
168 : !> \author fawzi (copying from the routine of Matthias Krack in functionals.F)
169 : !> \note
170 : !> Literature:- F. A. Hamprecht, A. J. Cohen, D. J. Tozer, and N. C. Handy,
171 : !> J. Chem. Phys. 109, 6264 (1998) -> HCTH/93
172 : !> - A. D. Boese, N. L. Doltsinis, N. C. Handy, and M. Sprik,
173 : !> J. Chem. Phys. 112, 1670 (2000) -> HCTH/120 and HCTH/147
174 : !> - A. D. Boese and N. C. Handy,
175 : !> J. Chem. Phys. 114, 5497 (2001) -> HCTH/407
176 : !> - J. P. Perdew and Y. Wang,
177 : !> Phys. Rev. B 45, 13244 (1992) -> PW92
178 : ! **************************************************************************************************
179 565 : SUBROUTINE hcth_lda_calc(iparset, rho, norm_drho, e_0, e_rho, e_ndrho, &
180 : epsilon_rho, npoints)
181 : INTEGER, INTENT(IN) :: iparset
182 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho
183 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho
184 : REAL(kind=dp), INTENT(in) :: epsilon_rho
185 : INTEGER, INTENT(IN) :: npoints
186 :
187 : REAL(KIND=dp), DIMENSION(4), PARAMETER :: &
188 : beta0 = (/7.59570_dp, 3.58760_dp, 1.63820_dp, 0.49294_dp/), &
189 : beta1 = (/14.11890_dp, 6.19770_dp, 3.36620_dp, 0.62517_dp/)
190 : REAL(KIND=dp), PARAMETER :: a0 = 0.031091_dp, a1 = 0.015545_dp, alpha0 = 0.21370_dp, &
191 : alpha1 = 0.20548_dp, f13 = 1.0_dp/3.0_dp, f43 = 4.0_dp*f13, f83 = 8.0_dp*f13, &
192 : gamma_cab = 0.006_dp, gamma_css = 0.200_dp, gamma_xss = 0.004_dp
193 :
194 : INTEGER :: ii
195 : REAL(KIND=dp) :: cx_vwn_e, cx_vwn_v, dgcabddrho, dgcabdrho, dgcabds, dgcssddrho, dgcssdrho, &
196 : dgcssds, dgdrs, dgxssddrho, dgxssdrho, dgxssds, drho, drhos, drsdrho, ecab, ecss, exss, &
197 : g, gcab, gcss, gs2, gxss, my_rho, p, q, rho13, rho43, rhos, rhos13, rhos43, rs, rs12, &
198 : rsfac, s, s2, two13, u, vcab, vcss, vxss, x, y
199 : REAL(KIND=dp), DIMENSION(0:4) :: ccab, ccss, cxss
200 :
201 565 : cx_vwn_e = -0.75_dp*(3.0_dp/pi)**f13
202 565 : cx_vwn_v = f43*cx_vwn_e
203 565 : rsfac = (f43*pi)**(-f13)
204 565 : two13 = 2.0_dp**f13
205 :
206 : ! *** LSDA correlation parametrisation (PW92) ***
207 : ! *** GGA parametrisation (HCTH/iparset) ***
208 :
209 : ! *** Load the HCTH parameter set HCTH/iparset ***
210 :
211 565 : SELECT CASE (iparset)
212 : CASE (93)
213 0 : cxss(0) = 0.109320E+01_dp
214 0 : ccss(0) = 0.222601E+00_dp
215 0 : ccab(0) = 0.729974E+00_dp
216 0 : cxss(1) = -0.744056E+00_dp
217 0 : ccss(1) = -0.338622E-01_dp
218 0 : ccab(1) = 0.335287E+01_dp
219 0 : cxss(2) = 0.559920E+01_dp
220 0 : ccss(2) = -0.125170E-01_dp
221 0 : ccab(2) = -0.115430E+02_dp
222 0 : cxss(3) = -0.678549E+01_dp
223 0 : ccss(3) = -0.802496E+00_dp
224 0 : ccab(3) = 0.808564E+01_dp
225 0 : cxss(4) = 0.449357E+01_dp
226 0 : ccss(4) = 0.155396E+01_dp
227 0 : ccab(4) = -0.447857E+01_dp
228 : CASE (120)
229 565 : cxss(0) = 0.109163E+01_dp
230 565 : ccss(0) = 0.489508E+00_dp
231 565 : ccab(0) = 0.514730E+00_dp
232 565 : cxss(1) = -0.747215E+00_dp
233 565 : ccss(1) = -0.260699E+00_dp
234 565 : ccab(1) = 0.692982E+01_dp
235 565 : cxss(2) = 0.507833E+01_dp
236 565 : ccss(2) = 0.432917E+00_dp
237 565 : ccab(2) = -0.247073E+02_dp
238 565 : cxss(3) = -0.410746E+01_dp
239 565 : ccss(3) = -0.199247E+01_dp
240 565 : ccab(3) = 0.231098E+02_dp
241 565 : cxss(4) = 0.117173E+01_dp
242 565 : ccss(4) = 0.248531E+01_dp
243 565 : ccab(4) = -0.113234E+02_dp
244 : CASE (147)
245 0 : cxss(0) = 0.109025E+01_dp
246 0 : ccss(0) = 0.562576E+00_dp
247 0 : ccab(0) = 0.542352E+00_dp
248 0 : cxss(1) = -0.799194E+00_dp
249 0 : ccss(1) = 0.171436E-01_dp
250 0 : ccab(1) = 0.701464E+01_dp
251 0 : cxss(2) = 0.557212E+01_dp
252 0 : ccss(2) = -0.130636E+01_dp
253 0 : ccab(2) = -0.283822E+02_dp
254 0 : cxss(3) = -0.586760E+01_dp
255 0 : ccss(3) = 0.105747E+01_dp
256 0 : ccab(3) = 0.350329E+02_dp
257 0 : cxss(4) = 0.304544E+01_dp
258 0 : ccss(4) = 0.885429E+00_dp
259 0 : ccab(4) = -0.204284E+02_dp
260 : CASE (407)
261 0 : cxss(0) = 0.108184E+01_dp
262 0 : ccss(0) = 0.118777E+01_dp
263 0 : ccab(0) = 0.589076E+00_dp
264 0 : cxss(1) = -0.518339E+00_dp
265 0 : ccss(1) = -0.240292E+01_dp
266 0 : ccab(1) = 0.442374E+01_dp
267 0 : cxss(2) = 0.342562E+01_dp
268 0 : ccss(2) = 0.561741E+01_dp
269 0 : ccab(2) = -0.192218E+02_dp
270 0 : cxss(3) = -0.262901E+01_dp
271 0 : ccss(3) = -0.917923E+01_dp
272 0 : ccab(3) = 0.425721E+02_dp
273 0 : cxss(4) = 0.228855E+01_dp
274 0 : ccss(4) = 0.624798E+01_dp
275 0 : ccab(4) = -0.420052E+02_dp
276 : ! DMB all-in-one HLE16 and applying 5/4 scaling
277 : ! to all exchange terms and 1/2 to all correlation terms
278 : CASE (408)
279 0 : cxss(0) = 0.108184E+01_dp*1.25_dp
280 0 : ccss(0) = 0.118777E+01_dp*0.5_dp
281 0 : ccab(0) = 0.589076E+00_dp*0.5_dp
282 0 : cxss(1) = -0.518339E+00_dp*1.25_dp
283 0 : ccss(1) = -0.240292E+01_dp*0.5_dp
284 0 : ccab(1) = 0.442374E+01_dp*0.5_dp
285 0 : cxss(2) = 0.342562E+01_dp*1.25_dp
286 0 : ccss(2) = 0.561741E+01_dp*0.5_dp
287 0 : ccab(2) = -0.192218E+02_dp*0.5_dp
288 0 : cxss(3) = -0.262901E+01_dp*1.25_dp
289 0 : ccss(3) = -0.917923E+01_dp*0.5_dp
290 0 : ccab(3) = 0.425721E+02_dp*0.5_dp
291 0 : cxss(4) = 0.228855E+01_dp*1.25_dp
292 0 : ccss(4) = 0.624798E+01_dp*0.5_dp
293 0 : ccab(4) = -0.420052E+02_dp*0.5_dp
294 : CASE DEFAULT
295 565 : CPABORT("Invalid HCTH parameter set requested ("//cp_to_string(iparset)//")")
296 : END SELECT
297 :
298 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(rho,norm_drho,cxss,ccss,&
299 : !$OMP ccab,cx_vwn_e, cx_vwn_v, rsfac, two13,epsilon_rho,npoints, &
300 : !$OMP e_0,e_rho,e_ndrho)&
301 : !$OMP PRIVATE(ii, dgcabddrho, dgcabdrho, dgcabds, dgcssddrho, &
302 : !$OMP dgcssdrho, dgcssds, dgdrs, dgxssddrho, dgxssdrho, dgxssds,&
303 : !$OMP drhos, drsdrho, ecab, ecss, exss, g, gcab, gcss, gs2, &
304 : !$OMP gxss, p, q, rho13, rho43, rhos, rhos13, rhos43, rs, rs12,&
305 565 : !$OMP s, s2, u, vcab, vcss, vxss, x, y, my_rho, drho)
306 : DO ii = 1, npoints
307 : ! *** rho_sigma = rho/2 = rho_alpha = rho_beta (same for |nabla rho|) ***
308 :
309 : IF (rho(ii) > epsilon_rho) THEN
310 : my_rho = MAX(rho(ii), epsilon_rho)
311 : drho = norm_drho(ii)
312 : rhos = 0.5_dp*my_rho
313 : drhos = 0.5_dp*drho
314 :
315 : rhos13 = rhos**f13
316 : rhos43 = rhos13*rhos
317 :
318 : rho13 = two13*rhos13
319 : rho43 = rho13*my_rho
320 :
321 : ! *** LSDA exchange part (VWN) ***
322 :
323 : exss = cx_vwn_e*rho43
324 : vxss = cx_vwn_v*rho13
325 :
326 : ! *** LSDA correlation part (PW92) ***
327 :
328 : ! *** G(rho_sigma,0) => spin polarisation zeta = 1 ***
329 :
330 : rs = rsfac/rhos13
331 : rs12 = SQRT(rs)
332 : q = 2.0_dp*a1*(beta1(1) + (beta1(2) + (beta1(3) + &
333 : beta1(4)*rs12)*rs12)*rs12)*rs12
334 : p = 1.0_dp + 1.0_dp/q
335 : x = -2.0_dp*a1*(1.0_dp + alpha1*rs)
336 : y = LOG(p)
337 : g = x*y
338 : dgdrs = -2.0_dp*a1*alpha1*y - &
339 : x*a1*(beta1(1)/rs12 + 2.0_dp*beta1(2) + &
340 : 3.0_dp*beta1(3)*rs12 + 4.0_dp*beta1(4)*rs)/(p*q*q)
341 : drsdrho = -f13*rs/my_rho
342 : ecss = my_rho*g
343 : vcss = g + my_rho*dgdrs*drsdrho
344 :
345 : ! *** G(rho_alpha,rho_beta) => spin polarisation zeta = 0 ***
346 :
347 : rs = rsfac/rho13
348 : rs12 = SQRT(rs)
349 : q = 2.0_dp*a0*(beta0(1) + (beta0(2) + (beta0(3) + &
350 : beta0(4)*rs12)*rs12)*rs12)*rs12
351 : p = 1.0_dp + 1.0_dp/q
352 : x = -2.0_dp*a0*(1.0_dp + alpha0*rs)
353 : y = LOG(p)
354 : g = x*y
355 : dgdrs = -2.0_dp*a0*alpha0*y - &
356 : x*a0*(beta0(1)/rs12 + 2.0_dp*beta0(2) + &
357 : 3.0_dp*beta0(3)*rs12 + 4.0_dp*beta0(4)*rs)/(p*q*q)
358 : drsdrho = -f13*rs/my_rho
359 : ecab = my_rho*g - ecss
360 : vcab = g + my_rho*dgdrs*drsdrho - vcss
361 :
362 : ! *** GGA part (HCTH) ***
363 :
364 : s = drhos/rhos43
365 : s2 = s*s
366 : x = -f83/my_rho
367 : y = 2.0_dp/(drho*drho)
368 :
369 : ! *** g_x(rho_sigma,rho_sigma) ***
370 :
371 : gs2 = gamma_xss*s2
372 : q = 1.0_dp/(1.0_dp + gs2)
373 : u = gs2*q
374 : gxss = cxss(0) + (cxss(1) + (cxss(2) + (cxss(3) + cxss(4)*u)*u)*u)*u
375 : dgxssds = q*(cxss(1) + (2.0_dp*cxss(2) + (3.0_dp*cxss(3) + &
376 : 4.0_dp*cxss(4)*u)*u)*u)*u
377 : dgxssdrho = x*dgxssds
378 : dgxssddrho = y*dgxssds
379 :
380 : ! *** g_c(rho_sigma,rho_sigma) ***
381 :
382 : gs2 = gamma_css*s2
383 : q = 1.0_dp/(1.0_dp + gs2)
384 : u = gs2*q
385 : gcss = ccss(0) + (ccss(1) + (ccss(2) + (ccss(3) + ccss(4)*u)*u)*u)*u
386 : dgcssds = q*(ccss(1) + (2.0_dp*ccss(2) + (3.0_dp*ccss(3) + &
387 : 4.0_dp*ccss(4)*u)*u)*u)*u
388 : dgcssdrho = x*dgcssds
389 : dgcssddrho = y*dgcssds
390 :
391 : ! *** g_c(rho_alpha,rho_beta) ***
392 :
393 : gs2 = gamma_cab*s2
394 : q = 1.0_dp/(1.0_dp + gs2)
395 : u = gs2*q
396 : gcab = ccab(0) + (ccab(1) + (ccab(2) + (ccab(3) + ccab(4)*u)*u)*u)*u
397 : dgcabds = q*(ccab(1) + (2.0_dp*ccab(2) + (3.0_dp*ccab(3) + &
398 : 4.0_dp*ccab(4)*u)*u)*u)*u
399 : dgcabdrho = x*dgcabds
400 : dgcabddrho = y*dgcabds
401 :
402 : ! *** Finally collect all contributions ***
403 :
404 : e_0(ii) = e_0(ii) + exss*gxss + ecss*gcss + ecab*gcab
405 : e_rho(ii) = e_rho(ii) + vxss*gxss + exss*dgxssdrho + &
406 : vcss*gcss + ecss*dgcssdrho + &
407 : vcab*gcab + ecab*dgcabdrho
408 : e_ndrho(ii) = e_ndrho(ii) + (exss*dgxssddrho + ecss*dgcssddrho + ecab*dgcabddrho)*drho
409 : END IF
410 : END DO
411 :
412 565 : END SUBROUTINE hcth_lda_calc
413 :
414 : END MODULE xc_hcth
|
