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 optx
10 : !> \note
11 : !> will need proper testing / review
12 : !> \author Joost VandeVondele [03.2004]
13 : ! **************************************************************************************************
14 : MODULE xc_optx
15 : USE cp_array_utils, ONLY: cp_3d_r_cp_type
16 : USE input_section_types, ONLY: section_vals_type,&
17 : section_vals_val_get
18 : USE kinds, ONLY: dp
19 : USE xc_derivative_desc, ONLY: deriv_norm_drho,&
20 : deriv_norm_drhoa,&
21 : deriv_norm_drhob,&
22 : deriv_rho,&
23 : deriv_rhoa,&
24 : deriv_rhob
25 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
26 : xc_dset_get_derivative
27 : USE xc_derivative_types, ONLY: xc_derivative_get,&
28 : xc_derivative_type
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 : PRIVATE
36 :
37 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_optx'
38 :
39 : PUBLIC :: optx_lda_info, optx_lda_eval, optx_lsd_info, optx_lsd_eval
40 : CONTAINS
41 :
42 : ! **************************************************************************************************
43 : !> \brief info about the optx functional
44 : !> \param reference string with the reference of the actual functional
45 : !> \param shortform string with the shortform of the functional name
46 : !> \param needs the components needed by this functional are set to
47 : !> true (does not set the unneeded components to false)
48 : !> \param max_deriv implemented derivative of the xc functional
49 : !> \author Joost
50 : ! **************************************************************************************************
51 435 : SUBROUTINE optx_lda_info(reference, shortform, needs, max_deriv)
52 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
53 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
54 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
55 :
56 435 : IF (PRESENT(reference)) THEN
57 3 : reference = "OPTX, Handy NC and Cohen AJ, JCP 116, p. 5411 (2002) (LDA)"
58 : END IF
59 435 : IF (PRESENT(shortform)) THEN
60 3 : shortform = "OPTX exchange (LDA)"
61 : END IF
62 435 : IF (PRESENT(needs)) THEN
63 432 : needs%rho = .TRUE.
64 432 : needs%norm_drho = .TRUE.
65 : END IF
66 435 : IF (PRESENT(max_deriv)) max_deriv = 1
67 435 : END SUBROUTINE optx_lda_info
68 :
69 : ! **************************************************************************************************
70 : !> \brief info about the optx functional (LSD)
71 : !> \param reference string with the reference of the actual functional
72 : !> \param shortform string with the shortform of the functional name
73 : !> \param needs the components needed by this functional are set to
74 : !> true (does not set the unneeded components to false)
75 : !> \param max_deriv implemented derivative of the xc functional
76 : !> \author Joost
77 : ! **************************************************************************************************
78 481 : SUBROUTINE optx_lsd_info(reference, shortform, needs, max_deriv)
79 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
80 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
81 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
82 :
83 481 : IF (PRESENT(reference)) THEN
84 1 : reference = "OPTX, Handy NC and Cohen AJ, JCP 116, p. 5411 (2002), (LSD) "
85 : END IF
86 481 : IF (PRESENT(shortform)) THEN
87 1 : shortform = "OPTX exchange (LSD)"
88 : END IF
89 481 : IF (PRESENT(needs)) THEN
90 480 : needs%rho_spin = .TRUE.
91 480 : needs%norm_drho_spin = .TRUE.
92 : END IF
93 481 : IF (PRESENT(max_deriv)) max_deriv = 1
94 481 : END SUBROUTINE optx_lsd_info
95 :
96 : ! **************************************************************************************************
97 : !> \brief evaluates the optx functional for lda
98 : !> \param rho_set the density where you want to evaluate the functional
99 : !> \param deriv_set place where to store the functional derivatives (they are
100 : !> added to the derivatives)
101 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
102 : !> \param optx_params input parameter (scaling)
103 : !> \par History
104 : !> 01.2007 added scaling [Manuel Guidon]
105 : !> \author Joost
106 : ! **************************************************************************************************
107 2820 : SUBROUTINE optx_lda_eval(rho_set, deriv_set, grad_deriv, optx_params)
108 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
109 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
110 : INTEGER, INTENT(in) :: grad_deriv
111 : TYPE(section_vals_type), POINTER :: optx_params
112 :
113 : INTEGER :: npoints
114 : INTEGER, DIMENSION(2, 3) :: bo
115 : REAL(kind=dp) :: a1, a2, epsilon_drho, epsilon_rho, gam, &
116 : sx
117 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
118 564 : POINTER :: e_0, e_ndrho, e_rho, norm_drho, rho
119 : TYPE(xc_derivative_type), POINTER :: deriv
120 :
121 564 : NULLIFY (e_0, e_ndrho, e_rho, norm_drho, rho)
122 :
123 564 : CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
124 564 : CALL section_vals_val_get(optx_params, "a1", r_val=a1)
125 564 : CALL section_vals_val_get(optx_params, "a2", r_val=a2)
126 564 : CALL section_vals_val_get(optx_params, "gamma", r_val=gam)
127 :
128 : CALL xc_rho_set_get(rho_set, rho=rho, &
129 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
130 564 : drho_cutoff=epsilon_drho)
131 564 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
132 :
133 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
134 564 : allocate_deriv=.TRUE.)
135 564 : CALL xc_derivative_get(deriv, deriv_data=e_0)
136 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
137 564 : allocate_deriv=.TRUE.)
138 564 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
139 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
140 564 : allocate_deriv=.TRUE.)
141 564 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
142 564 : IF (grad_deriv > 1) THEN
143 0 : CPABORT("derivatives bigger than 1 not implemented")
144 : END IF
145 :
146 : CALL optx_lda_calc(rho=rho, norm_drho=norm_drho, &
147 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
148 : npoints=npoints, epsilon_rho=epsilon_rho, &
149 : epsilon_drho=epsilon_drho, sx=sx, &
150 564 : a1=a1, a2=a2, gam=gam)
151 564 : END SUBROUTINE optx_lda_eval
152 :
153 : ! **************************************************************************************************
154 : !> \brief evaluates the optx functional for lsd
155 : !> \param rho_set the density where you want to evaluate the functional
156 : !> \param deriv_set place where to store the functional derivatives (they are
157 : !> added to the derivatives)
158 : !> \param grad_deriv maximum degree of the derivative that should be evaluated
159 : !> \param optx_params input parameter (scaling)
160 : !> \par History
161 : !> 01.2007 added scaling [Manuel Guidon]
162 : !> \author Joost
163 : ! **************************************************************************************************
164 2856 : SUBROUTINE optx_lsd_eval(rho_set, deriv_set, grad_deriv, optx_params)
165 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
166 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
167 : INTEGER, INTENT(in) :: grad_deriv
168 : TYPE(section_vals_type), POINTER :: optx_params
169 :
170 : INTEGER :: ispin, npoints
171 : INTEGER, DIMENSION(2, 3) :: bo
172 : REAL(kind=dp) :: a1, a2, epsilon_drho, epsilon_rho, gam, &
173 : sx
174 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
175 476 : POINTER :: e_0
176 5712 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_rho, ndrho, rho
177 : TYPE(xc_derivative_type), POINTER :: deriv
178 :
179 476 : NULLIFY (e_0)
180 1428 : DO ispin = 1, 2
181 952 : NULLIFY (e_rho(ispin)%array)
182 952 : NULLIFY (e_ndrho(ispin)%array)
183 952 : NULLIFY (rho(ispin)%array)
184 1428 : NULLIFY (ndrho(ispin)%array)
185 : END DO
186 :
187 476 : CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
188 476 : CALL section_vals_val_get(optx_params, "a1", r_val=a1)
189 476 : CALL section_vals_val_get(optx_params, "a2", r_val=a2)
190 476 : CALL section_vals_val_get(optx_params, "gamma", r_val=gam)
191 :
192 : CALL xc_rho_set_get(rho_set, rhoa=rho(1)%array, rhob=rho(2)%array, &
193 : norm_drhoa=ndrho(1)%array, &
194 : norm_drhob=ndrho(2)%array, rho_cutoff=epsilon_rho, &
195 476 : drho_cutoff=epsilon_drho, local_bounds=bo)
196 476 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
197 :
198 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
199 476 : allocate_deriv=.TRUE.)
200 476 : CALL xc_derivative_get(deriv, deriv_data=e_0)
201 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
202 476 : allocate_deriv=.TRUE.)
203 476 : CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
204 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
205 476 : allocate_deriv=.TRUE.)
206 476 : CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
207 :
208 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
209 476 : allocate_deriv=.TRUE.)
210 476 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
211 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
212 476 : allocate_deriv=.TRUE.)
213 476 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
214 :
215 476 : IF (grad_deriv > 1) THEN
216 0 : CPABORT("derivatives bigger than 1 not implemented")
217 : END IF
218 1428 : DO ispin = 1, 2
219 : CALL optx_lsd_calc(rho=rho(ispin)%array, norm_drho=ndrho(ispin)%array, &
220 : e_0=e_0, e_rho=e_rho(ispin)%array, e_ndrho=e_ndrho(ispin)%array, &
221 : npoints=npoints, epsilon_rho=epsilon_rho, &
222 : epsilon_drho=epsilon_drho, sx=sx, &
223 1428 : a1=a1, a2=a2, gam=gam)
224 : END DO
225 476 : END SUBROUTINE optx_lsd_eval
226 :
227 : ! **************************************************************************************************
228 : !> \brief optx exchange functional
229 : !> \param rho the full density
230 : !> \param norm_drho the norm of the gradient of the full density
231 : !> \param e_0 the value of the functional in that point
232 : !> \param e_rho the derivative of the functional wrt. rho
233 : !> \param e_ndrho the derivative of the functional wrt. norm_drho
234 : !> \param epsilon_rho the cutoff on rho
235 : !> \param epsilon_drho ...
236 : !> \param npoints ...
237 : !> \param sx scaling-parameter for exchange
238 : !> \param a1 a1 coefficient of the OPTX functional
239 : !> \param a2 a2 coefficient of the OPTX functional
240 : !> \param gam gamma coefficient of the OPTX functional
241 : !> \par History
242 : !> 01.2007 added scaling [Manuel Guidon]
243 : !> \author Joost VandeVondele
244 : ! **************************************************************************************************
245 564 : SUBROUTINE optx_lda_calc(rho, norm_drho, e_0, e_rho, e_ndrho, &
246 : epsilon_rho, epsilon_drho, npoints, sx, &
247 : a1, a2, gam)
248 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho
249 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho
250 : REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho
251 : INTEGER, INTENT(in) :: npoints
252 : REAL(kind=dp), INTENT(in) :: sx, a1, a2, gam
253 :
254 : REAL(KIND=dp), PARAMETER :: cx = 0.930525736349100_dp, &
255 : o43 = 4.0_dp/3.0_dp
256 :
257 : INTEGER :: ii
258 : REAL(KIND=dp) :: denom, ex, gamxsxs, myndrho, myrho, &
259 : rho43, tmp, xs
260 :
261 : !$OMP PARALLEL DO DEFAULT (NONE) &
262 : !$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
263 : !$OMP SHARED(epsilon_rho, epsilon_drho, sx, npoints) &
264 : !$OMP SHARED(a1, a2, gam) &
265 : !$OMP PRIVATE(ii, myrho, myndrho, rho43, xs, gamxsxs) &
266 564 : !$OMP PRIVATE(denom, ex, tmp)
267 :
268 : DO ii = 1, npoints
269 : ! we get the full density and need spin parts -> 0.5
270 : myrho = 0.5_dp*rho(ii)
271 : myndrho = 0.5_dp*MAX(norm_drho(ii), epsilon_drho)
272 : IF (myrho > 0.5_dp*epsilon_rho) THEN
273 : rho43 = myrho**o43
274 : xs = (myndrho/rho43)
275 : gamxsxs = gam*xs*xs
276 : denom = 1.0_dp/(1.0_dp + gamxsxs)
277 : ex = rho43*(a1*cx + a2*(gamxsxs*denom)**2)
278 : ! 2.0 for both spins
279 : e_0(ii) = e_0(ii) - (2.0_dp*ex)*sx
280 : tmp = rho43*2.0_dp*a2*gamxsxs*denom**2*(1.0_dp - gamxsxs*denom)
281 : ! derive e_0 wrt to rho (full) and ndrho (also full)
282 : e_rho(ii) = e_rho(ii) - ((o43*ex + tmp*gamxsxs*(-2.0_dp*o43))/myrho)*sx
283 : e_ndrho(ii) = e_ndrho(ii) - ((tmp*gam*2.0_dp*myndrho/rho43**2))*sx
284 : END IF
285 : END DO
286 :
287 : !$OMP END PARALLEL DO
288 :
289 564 : END SUBROUTINE optx_lda_calc
290 :
291 : ! **************************************************************************************************
292 : !> \brief optx exchange functional
293 : !> \param rho the *spin* density
294 : !> \param norm_drho the norm of the gradient of the *spin* density
295 : !> \param e_0 the value of the functional in that point
296 : !> \param e_rho the derivative of the functional wrt. rho
297 : !> \param e_ndrho the derivative of the functional wrt. norm_drho
298 : !> \param epsilon_rho the cutoff on rho
299 : !> \param epsilon_drho ...
300 : !> \param npoints ...
301 : !> \param sx scaling parameter for exchange
302 : !> \param a1 a1 coefficient of the OPTX functional
303 : !> \param a2 a2 coefficient of the OPTX functional
304 : !> \param gam gamma coefficient of the OPTX functional
305 : !> \par History
306 : !> 01.2007 added scaling [Manuel Guidon]
307 : !> \author Joost VandeVondele
308 : ! **************************************************************************************************
309 952 : SUBROUTINE optx_lsd_calc(rho, norm_drho, e_0, e_rho, e_ndrho, &
310 : epsilon_rho, epsilon_drho, npoints, sx, &
311 : a1, a2, gam)
312 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho
313 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho
314 : REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho
315 : INTEGER, INTENT(in) :: npoints
316 : REAL(kind=dp), INTENT(in) :: sx, a1, a2, gam
317 :
318 : REAL(KIND=dp), PARAMETER :: cx = 0.930525736349100_dp, &
319 : o43 = 4.0_dp/3.0_dp
320 :
321 : INTEGER :: ii
322 : REAL(KIND=dp) :: denom, ex, gamxsxs, myndrho, myrho, &
323 : rho43, tmp, xs
324 :
325 : !$OMP PARALLEL DO DEFAULT(NONE) &
326 : !$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
327 : !$OMP SHARED(epsilon_rho, epsilon_drho, npoints, sx) &
328 : !$OMP SHARED(a1, a2, gam) &
329 : !$OMP PRIVATE(ii, denom, ex, gamxsxs, myndrho, myrho) &
330 952 : !$OMP PRIVATE(rho43, tmp, xs)
331 :
332 : DO ii = 1, npoints
333 : ! we do have the spin density already
334 : myrho = rho(ii)
335 : myndrho = MAX(norm_drho(ii), epsilon_drho)
336 : IF (myrho > epsilon_rho) THEN
337 : rho43 = myrho**o43
338 : xs = (myndrho/rho43)
339 : gamxsxs = gam*xs*xs
340 : denom = 1.0_dp/(1.0_dp + gamxsxs)
341 : ex = rho43*(a1*cx + a2*(gamxsxs*denom)**2)
342 : ! for a single spin
343 : e_0(ii) = e_0(ii) - ex*sx
344 : tmp = rho43*2.0_dp*a2*gamxsxs*denom**2*(1.0_dp - gamxsxs*denom)
345 : ! derive e_0 wrt to rho and ndrho
346 : e_rho(ii) = e_rho(ii) - ((o43*ex + tmp*gamxsxs*(-2.0_dp*o43))/myrho)*sx
347 : e_ndrho(ii) = e_ndrho(ii) - ((tmp*gam*2.0_dp*myndrho/rho43**2))*sx
348 : END IF
349 : END DO
350 :
351 : !$OMP END PARALLEL DO
352 :
353 952 : END SUBROUTINE optx_lsd_calc
354 :
355 : END MODULE xc_optx
|
