Line data Source code
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 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 maximum degree of the derivative that should be evaluated
113 : !> \author fawzi
114 : ! **************************************************************************************************
115 565 : SUBROUTINE hcth_lda_eval(iparset, rho_set, deriv_set, grad_deriv)
116 : INTEGER, INTENT(in) :: iparset
117 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
118 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
119 : INTEGER, INTENT(in) :: grad_deriv
120 :
121 : INTEGER :: npoints
122 : INTEGER, DIMENSION(2, 3) :: bo
123 : REAL(kind=dp) :: epsilon_rho
124 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
125 565 : POINTER :: e_0, e_ndrho, e_rho, norm_drho, rho
126 : TYPE(xc_derivative_type), POINTER :: deriv
127 :
128 565 : NULLIFY (e_0, e_ndrho, e_rho, norm_drho, rho)
129 :
130 : CALL xc_rho_set_get(rho_set, rho=rho, &
131 565 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho)
132 565 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
133 :
134 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
135 565 : allocate_deriv=.TRUE.)
136 565 : CALL xc_derivative_get(deriv, deriv_data=e_0)
137 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
138 565 : allocate_deriv=.TRUE.)
139 565 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
140 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
141 565 : allocate_deriv=.TRUE.)
142 565 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
143 565 : IF (grad_deriv > 1) THEN
144 0 : CPABORT("derivatives bigger than 1 not implemented")
145 : END IF
146 :
147 : CALL hcth_lda_calc(iparset=iparset, rho=rho, norm_drho=norm_drho, &
148 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
149 565 : npoints=npoints, epsilon_rho=epsilon_rho)
150 565 : END SUBROUTINE hcth_lda_eval
151 :
152 : ! **************************************************************************************************
153 : !> \brief Calculate the gradient-corrected xc energy and potential
154 : !> of Hamprecht, Cohen, Tozer, and Handy (HCTH) for a closed shell
155 : !> density.
156 : !> \param iparset the parameter set that should be used (93,120,147,407)
157 : !> \param rho the density
158 : !> \param norm_drho the norm of the gradient of the density
159 : !> \param e_0 the value of the functional in that point
160 : !> \param e_rho the derivative of the functional wrt. rho
161 : !> \param e_ndrho the derivative of the functional wrt. norm_drho
162 : !> \param epsilon_rho the cutoff on rho
163 : !> \param npoints ...
164 : !> \author fawzi (copying from the routine of Matthias Krack in functionals.F)
165 : !> \note
166 : !> Literature:- F. A. Hamprecht, A. J. Cohen, D. J. Tozer, and N. C. Handy,
167 : !> J. Chem. Phys. 109, 6264 (1998) -> HCTH/93
168 : !> - A. D. Boese, N. L. Doltsinis, N. C. Handy, and M. Sprik,
169 : !> J. Chem. Phys. 112, 1670 (2000) -> HCTH/120 and HCTH/147
170 : !> - A. D. Boese and N. C. Handy,
171 : !> J. Chem. Phys. 114, 5497 (2001) -> HCTH/407
172 : !> - J. P. Perdew and Y. Wang,
173 : !> Phys. Rev. B 45, 13244 (1992) -> PW92
174 : ! **************************************************************************************************
175 565 : SUBROUTINE hcth_lda_calc(iparset, rho, norm_drho, e_0, e_rho, e_ndrho, &
176 : epsilon_rho, npoints)
177 : INTEGER, INTENT(IN) :: iparset
178 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho
179 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho
180 : REAL(kind=dp), INTENT(in) :: epsilon_rho
181 : INTEGER, INTENT(IN) :: npoints
182 :
183 : REAL(KIND=dp), DIMENSION(4), PARAMETER :: &
184 : beta0 = [7.59570_dp, 3.58760_dp, 1.63820_dp, 0.49294_dp], &
185 : beta1 = [14.11890_dp, 6.19770_dp, 3.36620_dp, 0.62517_dp]
186 : REAL(KIND=dp), PARAMETER :: a0 = 0.031091_dp, a1 = 0.015545_dp, alpha0 = 0.21370_dp, &
187 : alpha1 = 0.20548_dp, f13 = 1.0_dp/3.0_dp, f43 = 4.0_dp*f13, f83 = 8.0_dp*f13, &
188 : gamma_cab = 0.006_dp, gamma_css = 0.200_dp, gamma_xss = 0.004_dp
189 :
190 : INTEGER :: ii
191 : REAL(KIND=dp) :: cx_vwn_e, cx_vwn_v, dgcabddrho, dgcabdrho, dgcabds, dgcssddrho, dgcssdrho, &
192 : dgcssds, dgdrs, dgxssddrho, dgxssdrho, dgxssds, drho, drhos, drsdrho, ecab, ecss, exss, &
193 : g, gcab, gcss, gs2, gxss, my_rho, p, q, rho13, rho43, rhos, rhos13, rhos43, rs, rs12, &
194 : rsfac, s, s2, two13, u, vcab, vcss, vxss, x, y
195 : REAL(KIND=dp), DIMENSION(0:4) :: ccab, ccss, cxss
196 :
197 565 : cx_vwn_e = -0.75_dp*(3.0_dp/pi)**f13
198 565 : cx_vwn_v = f43*cx_vwn_e
199 565 : rsfac = (f43*pi)**(-f13)
200 565 : two13 = 2.0_dp**f13
201 :
202 : ! *** LSDA correlation parametrisation (PW92) ***
203 : ! *** GGA parametrisation (HCTH/iparset) ***
204 :
205 : ! *** Load the HCTH parameter set HCTH/iparset ***
206 :
207 565 : SELECT CASE (iparset)
208 : CASE (93)
209 0 : cxss(0) = 0.109320E+01_dp
210 0 : ccss(0) = 0.222601E+00_dp
211 0 : ccab(0) = 0.729974E+00_dp
212 0 : cxss(1) = -0.744056E+00_dp
213 0 : ccss(1) = -0.338622E-01_dp
214 0 : ccab(1) = 0.335287E+01_dp
215 0 : cxss(2) = 0.559920E+01_dp
216 0 : ccss(2) = -0.125170E-01_dp
217 0 : ccab(2) = -0.115430E+02_dp
218 0 : cxss(3) = -0.678549E+01_dp
219 0 : ccss(3) = -0.802496E+00_dp
220 0 : ccab(3) = 0.808564E+01_dp
221 0 : cxss(4) = 0.449357E+01_dp
222 0 : ccss(4) = 0.155396E+01_dp
223 0 : ccab(4) = -0.447857E+01_dp
224 : CASE (120)
225 565 : cxss(0) = 0.109163E+01_dp
226 565 : ccss(0) = 0.489508E+00_dp
227 565 : ccab(0) = 0.514730E+00_dp
228 565 : cxss(1) = -0.747215E+00_dp
229 565 : ccss(1) = -0.260699E+00_dp
230 565 : ccab(1) = 0.692982E+01_dp
231 565 : cxss(2) = 0.507833E+01_dp
232 565 : ccss(2) = 0.432917E+00_dp
233 565 : ccab(2) = -0.247073E+02_dp
234 565 : cxss(3) = -0.410746E+01_dp
235 565 : ccss(3) = -0.199247E+01_dp
236 565 : ccab(3) = 0.231098E+02_dp
237 565 : cxss(4) = 0.117173E+01_dp
238 565 : ccss(4) = 0.248531E+01_dp
239 565 : ccab(4) = -0.113234E+02_dp
240 : CASE (147)
241 0 : cxss(0) = 0.109025E+01_dp
242 0 : ccss(0) = 0.562576E+00_dp
243 0 : ccab(0) = 0.542352E+00_dp
244 0 : cxss(1) = -0.799194E+00_dp
245 0 : ccss(1) = 0.171436E-01_dp
246 0 : ccab(1) = 0.701464E+01_dp
247 0 : cxss(2) = 0.557212E+01_dp
248 0 : ccss(2) = -0.130636E+01_dp
249 0 : ccab(2) = -0.283822E+02_dp
250 0 : cxss(3) = -0.586760E+01_dp
251 0 : ccss(3) = 0.105747E+01_dp
252 0 : ccab(3) = 0.350329E+02_dp
253 0 : cxss(4) = 0.304544E+01_dp
254 0 : ccss(4) = 0.885429E+00_dp
255 0 : ccab(4) = -0.204284E+02_dp
256 : CASE (407)
257 0 : cxss(0) = 0.108184E+01_dp
258 0 : ccss(0) = 0.118777E+01_dp
259 0 : ccab(0) = 0.589076E+00_dp
260 0 : cxss(1) = -0.518339E+00_dp
261 0 : ccss(1) = -0.240292E+01_dp
262 0 : ccab(1) = 0.442374E+01_dp
263 0 : cxss(2) = 0.342562E+01_dp
264 0 : ccss(2) = 0.561741E+01_dp
265 0 : ccab(2) = -0.192218E+02_dp
266 0 : cxss(3) = -0.262901E+01_dp
267 0 : ccss(3) = -0.917923E+01_dp
268 0 : ccab(3) = 0.425721E+02_dp
269 0 : cxss(4) = 0.228855E+01_dp
270 0 : ccss(4) = 0.624798E+01_dp
271 0 : ccab(4) = -0.420052E+02_dp
272 : ! DMB all-in-one HLE16 and applying 5/4 scaling
273 : ! to all exchange terms and 1/2 to all correlation terms
274 : CASE (408)
275 0 : cxss(0) = 0.108184E+01_dp*1.25_dp
276 0 : ccss(0) = 0.118777E+01_dp*0.5_dp
277 0 : ccab(0) = 0.589076E+00_dp*0.5_dp
278 0 : cxss(1) = -0.518339E+00_dp*1.25_dp
279 0 : ccss(1) = -0.240292E+01_dp*0.5_dp
280 0 : ccab(1) = 0.442374E+01_dp*0.5_dp
281 0 : cxss(2) = 0.342562E+01_dp*1.25_dp
282 0 : ccss(2) = 0.561741E+01_dp*0.5_dp
283 0 : ccab(2) = -0.192218E+02_dp*0.5_dp
284 0 : cxss(3) = -0.262901E+01_dp*1.25_dp
285 0 : ccss(3) = -0.917923E+01_dp*0.5_dp
286 0 : ccab(3) = 0.425721E+02_dp*0.5_dp
287 0 : cxss(4) = 0.228855E+01_dp*1.25_dp
288 0 : ccss(4) = 0.624798E+01_dp*0.5_dp
289 0 : ccab(4) = -0.420052E+02_dp*0.5_dp
290 : CASE DEFAULT
291 565 : CPABORT("Invalid HCTH parameter set requested ("//cp_to_string(iparset)//")")
292 : END SELECT
293 :
294 : !$OMP PARALLEL DO DEFAULT(NONE) SHARED(rho,norm_drho,cxss,ccss,&
295 : !$OMP ccab,cx_vwn_e, cx_vwn_v, rsfac, two13,epsilon_rho,npoints, &
296 : !$OMP e_0,e_rho,e_ndrho)&
297 : !$OMP PRIVATE(ii, dgcabddrho, dgcabdrho, dgcabds, dgcssddrho, &
298 : !$OMP dgcssdrho, dgcssds, dgdrs, dgxssddrho, dgxssdrho, dgxssds,&
299 : !$OMP drhos, drsdrho, ecab, ecss, exss, g, gcab, gcss, gs2, &
300 : !$OMP gxss, p, q, rho13, rho43, rhos, rhos13, rhos43, rs, rs12,&
301 565 : !$OMP s, s2, u, vcab, vcss, vxss, x, y, my_rho, drho)
302 : DO ii = 1, npoints
303 : ! *** rho_sigma = rho/2 = rho_alpha = rho_beta (same for |nabla rho|) ***
304 :
305 : IF (rho(ii) > epsilon_rho) THEN
306 : my_rho = MAX(rho(ii), epsilon_rho)
307 : drho = norm_drho(ii)
308 : rhos = 0.5_dp*my_rho
309 : drhos = 0.5_dp*drho
310 :
311 : rhos13 = rhos**f13
312 : rhos43 = rhos13*rhos
313 :
314 : rho13 = two13*rhos13
315 : rho43 = rho13*my_rho
316 :
317 : ! *** LSDA exchange part (VWN) ***
318 :
319 : exss = cx_vwn_e*rho43
320 : vxss = cx_vwn_v*rho13
321 :
322 : ! *** LSDA correlation part (PW92) ***
323 :
324 : ! *** G(rho_sigma,0) => spin polarisation zeta = 1 ***
325 :
326 : rs = rsfac/rhos13
327 : rs12 = SQRT(rs)
328 : q = 2.0_dp*a1*(beta1(1) + (beta1(2) + (beta1(3) + &
329 : beta1(4)*rs12)*rs12)*rs12)*rs12
330 : p = 1.0_dp + 1.0_dp/q
331 : x = -2.0_dp*a1*(1.0_dp + alpha1*rs)
332 : y = LOG(p)
333 : g = x*y
334 : dgdrs = -2.0_dp*a1*alpha1*y - &
335 : x*a1*(beta1(1)/rs12 + 2.0_dp*beta1(2) + &
336 : 3.0_dp*beta1(3)*rs12 + 4.0_dp*beta1(4)*rs)/(p*q*q)
337 : drsdrho = -f13*rs/my_rho
338 : ecss = my_rho*g
339 : vcss = g + my_rho*dgdrs*drsdrho
340 :
341 : ! *** G(rho_alpha,rho_beta) => spin polarisation zeta = 0 ***
342 :
343 : rs = rsfac/rho13
344 : rs12 = SQRT(rs)
345 : q = 2.0_dp*a0*(beta0(1) + (beta0(2) + (beta0(3) + &
346 : beta0(4)*rs12)*rs12)*rs12)*rs12
347 : p = 1.0_dp + 1.0_dp/q
348 : x = -2.0_dp*a0*(1.0_dp + alpha0*rs)
349 : y = LOG(p)
350 : g = x*y
351 : dgdrs = -2.0_dp*a0*alpha0*y - &
352 : x*a0*(beta0(1)/rs12 + 2.0_dp*beta0(2) + &
353 : 3.0_dp*beta0(3)*rs12 + 4.0_dp*beta0(4)*rs)/(p*q*q)
354 : drsdrho = -f13*rs/my_rho
355 : ecab = my_rho*g - ecss
356 : vcab = g + my_rho*dgdrs*drsdrho - vcss
357 :
358 : ! *** GGA part (HCTH) ***
359 :
360 : s = drhos/rhos43
361 : s2 = s*s
362 : x = -f83/my_rho
363 : y = 2.0_dp/(drho*drho)
364 :
365 : ! *** g_x(rho_sigma,rho_sigma) ***
366 :
367 : gs2 = gamma_xss*s2
368 : q = 1.0_dp/(1.0_dp + gs2)
369 : u = gs2*q
370 : gxss = cxss(0) + (cxss(1) + (cxss(2) + (cxss(3) + cxss(4)*u)*u)*u)*u
371 : dgxssds = q*(cxss(1) + (2.0_dp*cxss(2) + (3.0_dp*cxss(3) + &
372 : 4.0_dp*cxss(4)*u)*u)*u)*u
373 : dgxssdrho = x*dgxssds
374 : dgxssddrho = y*dgxssds
375 :
376 : ! *** g_c(rho_sigma,rho_sigma) ***
377 :
378 : gs2 = gamma_css*s2
379 : q = 1.0_dp/(1.0_dp + gs2)
380 : u = gs2*q
381 : gcss = ccss(0) + (ccss(1) + (ccss(2) + (ccss(3) + ccss(4)*u)*u)*u)*u
382 : dgcssds = q*(ccss(1) + (2.0_dp*ccss(2) + (3.0_dp*ccss(3) + &
383 : 4.0_dp*ccss(4)*u)*u)*u)*u
384 : dgcssdrho = x*dgcssds
385 : dgcssddrho = y*dgcssds
386 :
387 : ! *** g_c(rho_alpha,rho_beta) ***
388 :
389 : gs2 = gamma_cab*s2
390 : q = 1.0_dp/(1.0_dp + gs2)
391 : u = gs2*q
392 : gcab = ccab(0) + (ccab(1) + (ccab(2) + (ccab(3) + ccab(4)*u)*u)*u)*u
393 : dgcabds = q*(ccab(1) + (2.0_dp*ccab(2) + (3.0_dp*ccab(3) + &
394 : 4.0_dp*ccab(4)*u)*u)*u)*u
395 : dgcabdrho = x*dgcabds
396 : dgcabddrho = y*dgcabds
397 :
398 : ! *** Finally collect all contributions ***
399 :
400 : e_0(ii) = e_0(ii) + exss*gxss + ecss*gcss + ecab*gcab
401 : e_rho(ii) = e_rho(ii) + vxss*gxss + exss*dgxssdrho + &
402 : vcss*gcss + ecss*dgcssdrho + &
403 : vcab*gcab + ecab*dgcabdrho
404 : e_ndrho(ii) = e_ndrho(ii) + (exss*dgxssddrho + ecss*dgcssddrho + ecab*dgcabddrho)*drho
405 : END IF
406 : END DO
407 :
408 565 : END SUBROUTINE hcth_lda_calc
409 :
410 : END MODULE xc_hcth
|
