Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 degree of the derivative that should be evaluated,
119 : !> if positive all the derivatives up to the given degree are evaluated,
120 : !> if negative only the given degree is calculated
121 : !> \param params input parameters (scaling, cutoff_radius)
122 : !> \par History
123 : !> 12.2008 created [Manuel Guidon]
124 : !> \author Manuel Guidon
125 : ! **************************************************************************************************
126 88 : SUBROUTINE xlda_hole_t_c_lr_lda_eval(rho_set, deriv_set, order, params)
127 :
128 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
129 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
130 : INTEGER, INTENT(IN) :: order
131 : TYPE(section_vals_type), POINTER :: params
132 :
133 : CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lda_eval'
134 :
135 : INTEGER :: handle, npoints
136 : INTEGER, DIMENSION(2, 3) :: bo
137 : REAL(kind=dp) :: epsilon_rho, R, sx
138 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
139 22 : POINTER :: dummy, e_0, e_rho, rho
140 : TYPE(xc_derivative_type), POINTER :: deriv
141 :
142 22 : CALL timeset(routineN, handle)
143 :
144 22 : CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
145 22 : CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
146 :
147 : CALL xc_rho_set_get(rho_set, rho=rho, &
148 22 : local_bounds=bo, rho_cutoff=epsilon_rho)
149 22 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
150 :
151 22 : dummy => rho
152 :
153 22 : e_0 => dummy
154 22 : e_rho => dummy
155 :
156 22 : IF (order >= 0) THEN
157 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
158 22 : allocate_deriv=.TRUE.)
159 22 : CALL xc_derivative_get(deriv, deriv_data=e_0)
160 : END IF
161 22 : IF (order >= 1 .OR. order == -1) THEN
162 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
163 22 : allocate_deriv=.TRUE.)
164 22 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
165 : END IF
166 22 : IF (order > 1 .OR. order < -1) THEN
167 0 : CPABORT("derivatives bigger than 1 not implemented")
168 : END IF
169 :
170 22 : IF (R == 0.0_dp) THEN
171 0 : CPABORT("Cutoff_Radius 0.0 not implemented")
172 : END IF
173 : CALL xlda_hole_t_c_lr_lda_calc(npoints, order, rho=rho, &
174 : e_0=e_0, e_rho=e_rho, &
175 : epsilon_rho=epsilon_rho, &
176 22 : sx=sx, R=R)
177 :
178 22 : CALL timestop(handle)
179 :
180 22 : END SUBROUTINE xlda_hole_t_c_lr_lda_eval
181 :
182 : ! **************************************************************************************************
183 : !> \brief Call low level routine
184 : !> \param npoints ...
185 : !> \param order ...
186 : !> \param rho ...
187 : !> \param e_0 ...
188 : !> \param e_rho ...
189 : !> \param epsilon_rho ...
190 : !> \param sx ...
191 : !> \param R ...
192 : !> \par History
193 : !> 12.2008 created [Manuel Guidon]
194 : !> \author Manuel Guidon
195 : ! **************************************************************************************************
196 22 : SUBROUTINE xlda_hole_t_c_lr_lda_calc(npoints, order, rho, e_0, e_rho, &
197 : epsilon_rho, sx, R)
198 :
199 : INTEGER, INTENT(in) :: npoints, order
200 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
201 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R
202 :
203 : INTEGER :: ip
204 : REAL(dp) :: my_rho
205 :
206 : !$OMP PARALLEL DO DEFAULT(NONE) &
207 : !$OMP SHARED(npoints, rho, epsilon_rho, order, e_0, e_rho) &
208 : !$OMP SHARED(sx, r) &
209 22 : !$OMP PRIVATE(ip, my_rho)
210 :
211 : DO ip = 1, npoints
212 : my_rho = MAX(rho(ip), 0.0_dp)
213 : IF (my_rho > epsilon_rho) THEN
214 : CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_0(ip), e_rho(ip), &
215 : sx, R)
216 : END IF
217 : END DO
218 :
219 : !$OMP END PARALLEL DO
220 :
221 22 : END SUBROUTINE xlda_hole_t_c_lr_lda_calc
222 :
223 : ! **************************************************************************************************
224 : !> \brief low level routine
225 : !> \param order ...
226 : !> \param rho ...
227 : !> \param e_0 ...
228 : !> \param e_rho ...
229 : !> \param sx ...
230 : !> \param R ...
231 : !> \par History
232 : !> 12.2008 created [Manuel Guidon]
233 : !> \author Manuel Guidon
234 : ! **************************************************************************************************
235 17330946 : SUBROUTINE xlda_hole_t_c_lr_lda_calc_0(order, rho, e_0, e_rho, &
236 : sx, R)
237 : INTEGER, INTENT(IN) :: order
238 : REAL(KIND=dp), INTENT(IN) :: rho
239 : REAL(kind=dp), INTENT(INOUT) :: e_0, e_rho
240 : REAL(KIND=dp), INTENT(IN) :: sx, R
241 :
242 : REAL(KIND=dp) :: t1, t12, t14, t15, t19, t2, t22, t23, &
243 : t24, t25, t3, t32, t33, t36, t4, t41, &
244 : t46, t5, t6, t62, t64, t67, t68, t7, &
245 : t82, t86, t9, t91, t95
246 :
247 17330946 : IF (order >= 0) THEN
248 17330946 : t1 = rho**2
249 17330946 : t2 = t1*pi
250 17330946 : t3 = 3**(0.1e1_dp/0.3e1_dp)
251 17330946 : t4 = pi**2
252 17330946 : t5 = t4*rho
253 17330946 : t6 = t5**(0.1e1_dp/0.3e1_dp)
254 17330946 : t7 = t6**2
255 17330946 : t9 = t3/t7
256 17330946 : t12 = LOG(R*t3*t6)
257 17330946 : t14 = R**2
258 17330946 : t15 = t14**2
259 17330946 : t19 = 0.1e1_dp/D
260 17330946 : t22 = t3**2
261 17330946 : t23 = t22*t7
262 17330946 : t24 = D*t14*t23
263 17330946 : t25 = EXP(-t24)
264 17330946 : t32 = 9 + 4*A*t14*t23
265 17330946 : t33 = LOG(t32)
266 17330946 : t36 = D**2
267 17330946 : t41 = expint(1, t24)
268 17330946 : t46 = 0.1e1_dp/t36
269 17330946 : t62 = LOG(0.2e1_dp)
270 17330946 : t64 = LOG(A)
271 : t67 = A*t12 + 0.3e1_dp/0.2e1_dp*E*t15*t3*t6*t5*t19*t25 &
272 : - A*t33/0.2e1_dp + E/t36/D*t25 + A*t41/0.2e1_dp + E*t14 &
273 : *t22*t7*t46*t25 + B*t19*t25/0.2e1_dp + C*t46*t25/0.2e1_dp &
274 : + C*t14*t22*t7*t19*t25/0.2e1_dp + A*t62 + A*t64 &
275 17330946 : /0.2e1_dp
276 17330946 : t68 = t9*t67
277 17330946 : e_0 = e_0 + (0.2e1_dp/0.3e1_dp*t2*t68)*sx
278 : END IF
279 17330946 : IF (order >= 1 .OR. order == -1) THEN
280 9091948 : t82 = A/rho
281 9091948 : t86 = t4**2
282 9091948 : t91 = A**2
283 9091948 : t95 = 0.1e1_dp/t6*t4
284 : e_rho = e_rho + (0.4e1_dp/0.3e1_dp*rho*pi*t68 - 0.4e1_dp/0.9e1_dp*t1*t4*pi &
285 : *t3/t7/t5*t67 + 0.2e1_dp/0.3e1_dp*t2*t9*(t82/0.3e1_dp - &
286 : 0.3e1_dp*E*t15*t14*t86*rho*t25 - 0.4e1_dp/0.3e1_dp*t91*t14 &
287 : *t22*t95/t32 - t82*t25/0.3e1_dp - B*t14*t22*t95*t25 &
288 9091948 : /0.3e1_dp - C*t15*t3*t6*t4*t25))*sx
289 : END IF
290 :
291 17330946 : END SUBROUTINE xlda_hole_t_c_lr_lda_calc_0
292 :
293 : ! **************************************************************************************************
294 : !> \brief evaluates the truncated lsd exchange hole. Calls the lda routine and
295 : !> applies spin scaling relation
296 : !> \param rho_set the density where you want to evaluate the functional
297 : !> \param deriv_set place where to store the functional derivatives (they are
298 : !> added to the derivatives)
299 : !> \param order degree of the derivative that should be evaluated,
300 : !> if positive all the derivatives up to the given degree are evaluated,
301 : !> if negative only the given degree is calculated
302 : !> \param params input parameters (scaling, cutoff_radius)
303 : !> \par History
304 : !> 12.2008 created [Manuel Guidon]
305 : !> \author Manuel Guidon
306 : ! **************************************************************************************************
307 232 : SUBROUTINE xlda_hole_t_c_lr_lsd_eval(rho_set, deriv_set, order, params)
308 :
309 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
310 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
311 : INTEGER, INTENT(IN) :: order
312 : TYPE(section_vals_type), POINTER :: params
313 :
314 : CHARACTER(len=*), PARAMETER :: routineN = 'xlda_hole_t_c_lr_lsd_eval'
315 :
316 : INTEGER :: handle, npoints
317 : INTEGER, DIMENSION(2, 3) :: bo
318 : REAL(kind=dp) :: epsilon_rho, R, sx
319 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
320 58 : POINTER :: dummy, e_0, e_rhoa, e_rhob, rhoa, rhob
321 : TYPE(xc_derivative_type), POINTER :: deriv
322 :
323 58 : CALL timeset(routineN, handle)
324 :
325 58 : CALL section_vals_val_get(params, "SCALE_X", r_val=sx)
326 58 : CALL section_vals_val_get(params, "CUTOFF_RADIUS", r_val=R)
327 :
328 : CALL xc_rho_set_get(rho_set, rhoa=rhoa, rhob=rhob, &
329 58 : local_bounds=bo, rho_cutoff=epsilon_rho)
330 58 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
331 :
332 58 : dummy => rhoa
333 :
334 58 : e_0 => dummy
335 58 : e_rhoa => dummy
336 58 : e_rhob => dummy
337 :
338 58 : IF (order >= 0) THEN
339 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
340 58 : allocate_deriv=.TRUE.)
341 58 : CALL xc_derivative_get(deriv, deriv_data=e_0)
342 : END IF
343 58 : IF (order >= 1 .OR. order == -1) THEN
344 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
345 30 : allocate_deriv=.TRUE.)
346 30 : CALL xc_derivative_get(deriv, deriv_data=e_rhoa)
347 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
348 30 : allocate_deriv=.TRUE.)
349 30 : CALL xc_derivative_get(deriv, deriv_data=e_rhob)
350 : END IF
351 58 : IF (order > 1 .OR. order < -1) THEN
352 0 : CPABORT("derivatives bigger than 2 not implemented")
353 : END IF
354 58 : IF (R == 0.0_dp) THEN
355 0 : CPABORT("Cutoff_Radius 0.0 not implemented")
356 : END IF
357 :
358 : !$OMP PARALLEL DEFAULT(NONE) &
359 : !$OMP SHARED(npoints, order, rhoa, e_0, e_rhoa, epsilon_rho) &
360 58 : !$OMP SHARED(sx, r,rhob, e_rhob)
361 :
362 : CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhoa, &
363 : e_0=e_0, e_rho=e_rhoa, &
364 : epsilon_rho=epsilon_rho, &
365 : sx=sx, R=R)
366 :
367 : CALL xlda_hole_t_c_lr_lsd_calc(npoints, order, rho=rhob, &
368 : e_0=e_0, e_rho=e_rhob, &
369 : epsilon_rho=epsilon_rho, &
370 : sx=sx, R=R)
371 : !$OMP END PARALLEL
372 :
373 58 : CALL timestop(handle)
374 :
375 58 : END SUBROUTINE xlda_hole_t_c_lr_lsd_eval
376 :
377 : ! **************************************************************************************************
378 : !> \brief low level routine
379 : !> \param npoints ...
380 : !> \param order ...
381 : !> \param rho ...
382 : !> \param e_0 ...
383 : !> \param e_rho ...
384 : !> \param epsilon_rho ...
385 : !> \param sx ...
386 : !> \param R ...
387 : !> \par History
388 : !> 12.2008 created [Manuel Guidon]
389 : !> \author Manuel Guidon
390 : ! **************************************************************************************************
391 116 : SUBROUTINE xlda_hole_t_c_lr_lsd_calc(npoints, order, rho, e_0, e_rho, &
392 : epsilon_rho, sx, R)
393 :
394 : INTEGER, INTENT(in) :: npoints, order
395 : REAL(kind=dp), DIMENSION(1:npoints), INTENT(inout) :: rho, e_0, e_rho
396 : REAL(kind=dp), INTENT(in) :: epsilon_rho, sx, R
397 :
398 : INTEGER :: ip
399 : REAL(dp) :: e_tmp, my_rho
400 :
401 : !$OMP DO
402 :
403 : DO ip = 1, npoints
404 5285250 : my_rho = 2.0_dp*MAX(rho(ip), 0.0_dp)
405 5285250 : IF (my_rho > epsilon_rho) THEN
406 5285250 : e_tmp = 0.0_dp
407 : CALL xlda_hole_t_c_lr_lda_calc_0(order, my_rho, e_tmp, e_rho(ip), &
408 5285250 : sx, R)
409 5285250 : e_0(ip) = e_0(ip) + 0.5_dp*e_tmp
410 : END IF
411 : END DO
412 :
413 : !$OMP END DO
414 :
415 116 : END SUBROUTINE xlda_hole_t_c_lr_lsd_calc
416 : END MODULE xc_xlda_hole_t_c_lr
417 :
|
