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 Calculates the lda exchange hole in a truncated coulomb potential.
10 : !> Can be used as longrange correction for truncated hfx calculations
11 : !> \par History
12 : !> Manuel Guidon (12.2008) : created
13 : !> \author Manuel Guidon (06.2008)
14 : ! **************************************************************************************************
15 :
16 : MODULE xc_xlda_hole_t_c_lr
17 :
18 : USE input_section_types, ONLY: section_vals_type,&
19 : section_vals_val_get
20 : USE kinds, ONLY: dp
21 : USE mathconstants, ONLY: pi
22 : USE mathlib, ONLY: expint
23 : USE xc_derivative_desc, ONLY: deriv_rho,&
24 : deriv_rhoa,&
25 : deriv_rhob
26 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
27 : xc_dset_get_derivative
28 : USE xc_derivative_types, ONLY: xc_derivative_get,&
29 : xc_derivative_type
30 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
31 : USE xc_rho_set_types, ONLY: xc_rho_set_get,&
32 : xc_rho_set_type
33 : #include "../base/base_uses.f90"
34 :
35 : IMPLICIT NONE
36 :
37 : PRIVATE
38 :
39 : ! *** Global parameters ***
40 :
41 : PUBLIC :: xlda_hole_t_c_lr_lda_eval, xlda_hole_t_c_lr_lda_info, &
42 : xlda_hole_t_c_lr_lsd_eval, xlda_hole_t_c_lr_lsd_info, &
43 : xlda_hole_t_c_lr_lda_calc_0
44 :
45 : REAL(KIND=dp), PARAMETER :: A = 1.0161144_dp, &
46 : B = -0.37170836_dp, &
47 : C = -0.077215461_dp, &
48 : D = 0.57786348_dp, &
49 : E = -0.051955731_dp
50 :
51 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_xlda_hole_t_c_lr'
52 :
53 : CONTAINS
54 :
55 : ! **************************************************************************************************
56 : !> \brief returns various information on the functional
57 : !> \param reference string with the reference of the actual functional
58 : !> \param shortform string with the shortform of the functional name
59 : !> \param needs the components needed by this functional are set to
60 : !> true (does not set the unneeded components to false)
61 : !> \param max_deriv controls the number of derivatives
62 : !> \par History
63 : !> 12.2008 created [mguidon]
64 : !> \author mguidon
65 : ! **************************************************************************************************
66 27 : SUBROUTINE xlda_hole_t_c_lr_lda_info(reference, shortform, needs, max_deriv)
67 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
68 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
69 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
70 :
71 27 : IF (PRESENT(reference)) THEN
72 1 : reference = "{LDA version}"
73 : END IF
74 27 : IF (PRESENT(shortform)) THEN
75 1 : shortform = "{LDA}"
76 : END IF
77 27 : IF (PRESENT(needs)) THEN
78 26 : needs%rho = .TRUE.
79 : END IF
80 27 : IF (PRESENT(max_deriv)) max_deriv = 1
81 :
82 27 : END SUBROUTINE xlda_hole_t_c_lr_lda_info
83 :
84 : ! **************************************************************************************************
85 : !> \brief returns various information on the functional
86 : !> \param reference string with the reference of the actual functional
87 : !> \param shortform string with the shortform of the functional name
88 : !> \param needs the components needed by this functional are set to
89 : !> true (does not set the unneeded components to false)
90 : !> \param max_deriv controls the number of derivatives
91 : !> \par History
92 : !> 12.2008 created [mguidon]
93 : !> \author mguidon
94 : ! **************************************************************************************************
95 63 : SUBROUTINE xlda_hole_t_c_lr_lsd_info(reference, shortform, needs, max_deriv)
96 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
97 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
98 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
99 :
100 63 : IF (PRESENT(reference)) THEN
101 1 : reference = "{LSD version}"
102 : END IF
103 63 : IF (PRESENT(shortform)) THEN
104 1 : shortform = "{LSD}"
105 : END IF
106 63 : IF (PRESENT(needs)) THEN
107 62 : needs%rho_spin = .TRUE.
108 : END IF
109 63 : IF (PRESENT(max_deriv)) max_deriv = 1
110 :
111 63 : END SUBROUTINE xlda_hole_t_c_lr_lsd_info
112 :
113 : ! **************************************************************************************************
114 : !> \brief evaluates the truncated lda exchange hole
115 : !> \param rho_set the density where you want to evaluate the functional
116 : !> \param deriv_set place where to store the functional derivatives (they are
117 : !> added to the derivatives)
118 : !> \param order maximum degree of the derivative that should be evaluated
119 : !> \param params input parameters (scaling, cutoff_radius)
120 : !> \par History
121 : !> 12.2008 created [Manuel Guidon]
122 : !> \author Manuel Guidon
123 : ! **************************************************************************************************
124 88 : SUBROUTINE xlda_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
125 :
126 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
127 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
128 : INTEGER, INTENT(IN) :: order
129 : TYPE(section_vals_type), POINTER :: params
130 :
131 : CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lda_eval'
132 :
133 : INTEGER :: handle, npoints
134 : INTEGER, DIMENSION(2, 3) :: bo
135 : REAL(kind=dp) :: epsilon_rho, R, sx
136 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
137 22 : POINTER :: dummy, e_0, e_rho, rho
138 : TYPE(xc_derivative_type), POINTER :: deriv
139 :
140 22 : CALL timeset(routineN, handle)
141 :
142 22 : CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
143 22 : CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
144 :
145 : CALL xc_rho_set_get(rho_set, rho=rho, &
146 22 : local_bounds=bo, rho_cutoff=epsilon_rho)
147 22 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
148 :
149 22 : dummy => rho
150 :
151 22 : e_0 => dummy
152 22 : e_rho => dummy
153 :
154 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
155 22 : allocate_deriv=.TRUE.)
156 22 : CALL xc_derivative_get(deriv, deriv_data=e_0)
157 22 : IF (order >= 1) THEN
158 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
159 22 : allocate_deriv=.TRUE.)
160 22 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
161 : END IF
162 22 : IF (order > 1) THEN
163 0 : CPABORT("derivatives bigger than 1 not implemented")
164 : END IF
165 :
166 22 : IF (R == 0.0_dp) THEN
167 0 : CPABORT("Cutoff_Radius 0.0 not implemented")
168 : END IF
169 : CALL xlda_hole_t_c_lr_lda_calc(npoints, order, rho=rho, &
170 : e_0=e_0, e_rho=e_rho, &
171 : epsilon_rho=epsilon_rho, &
172 22 : sx=sx, R=R)
173 :
174 22 : CALL timestop(handle)
175 :
176 22 : END SUBROUTINE xlda_hole_t_c_lr_lda_eval
177 :
178 : ! **************************************************************************************************
179 : !> \brief Call low level routine
180 : !> \param npoints ...
181 : !> \param order ...
182 : !> \param rho ...
183 : !> \param e_0 ...
184 : !> \param e_rho ...
185 : !> \param epsilon_rho ...
186 : !> \param sx ...
187 : !> \param R ...
188 : !> \par History
189 : !> 12.2008 created [Manuel Guidon]
190 : !> \author Manuel Guidon
191 : ! **************************************************************************************************
192 22 : SUBROUTINE xlda_hole_t_c_lr_lda_calc(npoints, order, rho, e_0, e_rho, &
193 : epsilon_rho, sx, R)
194 :
195 : INTEGER, INTENT(in) :: npoints, order
196 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
197 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R
198 :
199 : INTEGER :: ip
200 : REAL(dp) :: my_rho
201 :
202 : !$OMP PARALLEL DO DEFAULT(NONE) &
203 : !$OMP SHARED(npoints, rho, epsilon_rho, order, e_0, e_rho) &
204 : !$OMP SHARED(sx, r) &
205 22 : !$OMP PRIVATE(ip, my_rho)
206 :
207 : DO ip = 1, npoints
208 : my_rho = MAX(rho(ip), 0.0_dp)
209 : IF (my_rho > epsilon_rho) THEN
210 : CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_0(ip), e_rho(ip), &
211 : sx, R)
212 : END IF
213 : END DO
214 :
215 : !$OMP END PARALLEL DO
216 :
217 22 : END SUBROUTINE xlda_hole_t_c_lr_lda_calc
218 :
219 : ! **************************************************************************************************
220 : !> \brief low level routine
221 : !> \param order ...
222 : !> \param rho ...
223 : !> \param e_0 ...
224 : !> \param e_rho ...
225 : !> \param sx ...
226 : !> \param R ...
227 : !> \par History
228 : !> 12.2008 created [Manuel Guidon]
229 : !> \author Manuel Guidon
230 : ! **************************************************************************************************
231 17330946 : SUBROUTINE xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho, &
232 : sx, R)
233 : INTEGER, INTENT(IN) :: order
234 : REAL(KIND=dp), INTENT(IN) :: rho
235 : REAL(kind=dp), INTENT(INOUT) :: e_0, e_rho
236 : REAL(KIND=dp), INTENT(IN) :: sx, R
237 :
238 : REAL(KIND=dp) :: t1, t12, t14, t15, t19, t2, t22, t23, &
239 : t24, t25, t3, t32, t33, t36, t4, t41, &
240 : t46, t5, t6, t62, t64, t67, t68, t7, &
241 : t82, t86, t9, t91, t95
242 :
243 17330946 : t1 = rho**2
244 17330946 : t2 = t1*pi
245 17330946 : t3 = 3**(0.1e1_dp/0.3e1_dp)
246 17330946 : t4 = pi**2
247 17330946 : t5 = t4*rho
248 17330946 : t6 = t5**(0.1e1_dp/0.3e1_dp)
249 17330946 : t7 = t6**2
250 17330946 : t9 = t3/t7
251 17330946 : t12 = LOG(R*t3*t6)
252 17330946 : t14 = R**2
253 17330946 : t15 = t14**2
254 17330946 : t19 = 0.1e1_dp/D
255 17330946 : t22 = t3**2
256 17330946 : t23 = t22*t7
257 17330946 : t24 = D*t14*t23
258 17330946 : t25 = EXP(-t24)
259 17330946 : t32 = 9 + 4*A*t14*t23
260 17330946 : t33 = LOG(t32)
261 17330946 : t36 = D**2
262 17330946 : t41 = expint(1, t24)
263 17330946 : t46 = 0.1e1_dp/t36
264 17330946 : t62 = LOG(0.2e1_dp)
265 17330946 : t64 = LOG(A)
266 : t67 = A*t12 + 0.3e1_dp/0.2e1_dp*E*t15*t3*t6*t5*t19*t25 &
267 : - A*t33/0.2e1_dp + E/t36/D*t25 + A*t41/0.2e1_dp + E*t14 &
268 : *t22*t7*t46*t25 + B*t19*t25/0.2e1_dp + C*t46*t25/0.2e1_dp &
269 : + C*t14*t22*t7*t19*t25/0.2e1_dp + A*t62 + A*t64 &
270 17330946 : /0.2e1_dp
271 17330946 : t68 = t9*t67
272 17330946 : e_0 = e_0 + (0.2e1_dp/0.3e1_dp*t2*t68)*sx
273 17330946 : IF (order >= 1) THEN
274 9091948 : t82 = A/rho
275 9091948 : t86 = t4**2
276 9091948 : t91 = A**2
277 9091948 : t95 = 0.1e1_dp/t6*t4
278 : e_rho = e_rho + (0.4e1_dp/0.3e1_dp*rho*pi*t68 - 0.4e1_dp/0.9e1_dp*t1*t4*pi &
279 : *t3/t7/t5*t67 + 0.2e1_dp/0.3e1_dp*t2*t9*(t82/0.3e1_dp - &
280 : 0.3e1_dp*E*t15*t14*t86*rho*t25 - 0.4e1_dp/0.3e1_dp*t91*t14 &
281 : *t22*t95/t32 - t82*t25/0.3e1_dp - B*t14*t22*t95*t25 &
282 9091948 : /0.3e1_dp - C*t15*t3*t6*t4*t25))*sx
283 : END IF
284 :
285 17330946 : END SUBROUTINE xlda_hole_t_c_lr_lda_calc_0
286 :
287 : ! **************************************************************************************************
288 : !> \brief evaluates the truncated lsd exchange hole. Calls the lda routine and
289 : !> applies spin scaling relation
290 : !> \param rho_set the density where you want to evaluate the functional
291 : !> \param deriv_set place where to store the functional derivatives (they are
292 : !> added to the derivatives)
293 : !> \param order maximum degree of the derivative that should be evaluated
294 : !> \param params input parameters (scaling, cutoff_radius)
295 : !> \par History
296 : !> 12.2008 created [Manuel Guidon]
297 : !> \author Manuel Guidon
298 : ! **************************************************************************************************
299 232 : SUBROUTINE xlda_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params)
300 :
301 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
302 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
303 : INTEGER, INTENT(IN) :: order
304 : TYPE(section_vals_type), POINTER :: params
305 :
306 : CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lsd_eval'
307 :
308 : INTEGER :: handle, npoints
309 : INTEGER, DIMENSION(2, 3) :: bo
310 : REAL(kind=dp) :: epsilon_rho, R, sx
311 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
312 58 : POINTER :: dummy, e_0, e_rhoa, e_rhob, rhoa, rhob
313 : TYPE(xc_derivative_type), POINTER :: deriv
314 :
315 58 : CALL timeset(routineN, handle)
316 :
317 58 : CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
318 58 : CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
319 :
320 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
321 58 : local_bounds=bo, rho_cutoff=epsilon_rho)
322 58 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
323 :
324 58 : dummy => rhoa
325 :
326 58 : e_0 => dummy
327 58 : e_rhoa => dummy
328 58 : e_rhob => dummy
329 :
330 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
331 58 : allocate_deriv=.TRUE.)
332 58 : CALL xc_derivative_get(deriv, deriv_data=e_0)
333 58 : IF (order >= 1) THEN
334 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
335 30 : allocate_deriv=.TRUE.)
336 30 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
337 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
338 30 : allocate_deriv=.TRUE.)
339 30 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
340 : END IF
341 58 : IF (order > 1) THEN
342 0 : CPABORT("derivatives bigger than 2 not implemented")
343 : END IF
344 58 : IF (R == 0.0_dp) THEN
345 0 : CPABORT("Cutoff_Radius 0.0 not implemented")
346 : END IF
347 :
348 : !$OMP PARALLEL DEFAULT(NONE) &
349 : !$OMP SHARED(npoints, order, rhoa, e_0, e_rhoa, epsilon_rho) &
350 58 : !$OMP SHARED(sx, r,rhob, e_rhob)
351 :
352 : CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhoa, &
353 : e_0=e_0, e_rho=e_rhoa, &
354 : epsilon_rho=epsilon_rho, &
355 : sx=sx, R=R)
356 :
357 : CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhob, &
358 : e_0=e_0, e_rho=e_rhob, &
359 : epsilon_rho=epsilon_rho, &
360 : sx=sx, R=R)
361 : !$OMP END PARALLEL
362 :
363 58 : CALL timestop(handle)
364 :
365 58 : END SUBROUTINE xlda_hole_t_c_lr_lsd_eval
366 :
367 : ! **************************************************************************************************
368 : !> \brief low level routine
369 : !> \param npoints ...
370 : !> \param order ...
371 : !> \param rho ...
372 : !> \param e_0 ...
373 : !> \param e_rho ...
374 : !> \param epsilon_rho ...
375 : !> \param sx ...
376 : !> \param R ...
377 : !> \par History
378 : !> 12.2008 created [Manuel Guidon]
379 : !> \author Manuel Guidon
380 : ! **************************************************************************************************
381 116 : SUBROUTINE xlda_hole_t_c_lr_lsd_calc(npoints, order, rho, e_0, e_rho, &
382 : epsilon_rho, sx, R)
383 :
384 : INTEGER, INTENT(in) :: npoints, order
385 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
386 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R
387 :
388 : INTEGER :: ip
389 : REAL(dp) :: e_tmp, my_rho
390 :
391 : !$OMP DO
392 :
393 : DO ip = 1, npoints
394 5285250 : my_rho = 2.0_dp*MAX(rho(ip), 0.0_dp)
395 5285250 : IF (my_rho > epsilon_rho) THEN
396 5285250 : e_tmp = 0.0_dp
397 : CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_tmp, e_rho(ip), &
398 5285250 : sx, R)
399 5285250 : e_0(ip) = e_0(ip) + 0.5_dp*e_tmp
400 : END IF
401 : END DO
402 :
403 : !$OMP END DO
404 :
405 116 : END SUBROUTINE xlda_hole_t_c_lr_lsd_calc
406 : END MODULE xc_xlda_hole_t_c_lr
407 :
|
