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 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 degree of the derivative that should be evaluated,
102 : !> if positive all the derivatives up to the given degree are evaluated,
103 : !> if negative only the given degree is calculated
104 : !> \param optx_params input parameter (scaling)
105 : !> \par History
106 : !> 01.2007 added scaling [Manuel Guidon]
107 : !> \author Joost
108 : ! **************************************************************************************************
109 2820 : SUBROUTINE optx_lda_eval(rho_set, deriv_set, grad_deriv, optx_params)
110 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
111 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
112 : INTEGER, INTENT(in) :: grad_deriv
113 : TYPE(section_vals_type), POINTER :: optx_params
114 :
115 : INTEGER :: npoints
116 : INTEGER, DIMENSION(2, 3) :: bo
117 : REAL(kind=dp) :: a1, a2, epsilon_drho, epsilon_rho, gam, &
118 : sx
119 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
120 564 : POINTER :: e_0, e_ndrho, e_rho, norm_drho, rho
121 : TYPE(xc_derivative_type), POINTER :: deriv
122 :
123 564 : NULLIFY (e_0, e_ndrho, e_rho, norm_drho, rho)
124 :
125 564 : CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
126 564 : CALL section_vals_val_get(optx_params, "a1", r_val=a1)
127 564 : CALL section_vals_val_get(optx_params, "a2", r_val=a2)
128 564 : CALL section_vals_val_get(optx_params, "gamma", r_val=gam)
129 :
130 : CALL xc_rho_set_get(rho_set, rho=rho, &
131 : norm_drho=norm_drho, local_bounds=bo, rho_cutoff=epsilon_rho, &
132 564 : drho_cutoff=epsilon_drho)
133 564 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
134 :
135 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
136 564 : allocate_deriv=.TRUE.)
137 564 : CALL xc_derivative_get(deriv, deriv_data=e_0)
138 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
139 564 : allocate_deriv=.TRUE.)
140 564 : CALL xc_derivative_get(deriv, deriv_data=e_rho)
141 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
142 564 : allocate_deriv=.TRUE.)
143 564 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
144 564 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
145 0 : CPABORT("derivatives bigger than 1 not implemented")
146 : END IF
147 :
148 : CALL optx_lda_calc(rho=rho, norm_drho=norm_drho, &
149 : e_0=e_0, e_rho=e_rho, e_ndrho=e_ndrho, &
150 : npoints=npoints, epsilon_rho=epsilon_rho, &
151 : epsilon_drho=epsilon_drho, sx=sx, &
152 564 : a1=a1, a2=a2, gam=gam)
153 564 : END SUBROUTINE optx_lda_eval
154 :
155 : ! **************************************************************************************************
156 : !> \brief evaluates the optx functional for lsd
157 : !> \param rho_set the density where you want to evaluate the functional
158 : !> \param deriv_set place where to store the functional derivatives (they are
159 : !> added to the derivatives)
160 : !> \param grad_deriv degree of the derivative that should be evaluated,
161 : !> if positive all the derivatives up to the given degree are evaluated,
162 : !> if negative only the given degree is calculated
163 : !> \param optx_params input parameter (scaling)
164 : !> \par History
165 : !> 01.2007 added scaling [Manuel Guidon]
166 : !> \author Joost
167 : ! **************************************************************************************************
168 2856 : SUBROUTINE optx_lsd_eval(rho_set, deriv_set, grad_deriv, optx_params)
169 : TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
170 : TYPE(xc_derivative_set_type), INTENT(IN) :: deriv_set
171 : INTEGER, INTENT(in) :: grad_deriv
172 : TYPE(section_vals_type), POINTER :: optx_params
173 :
174 : INTEGER :: ispin, npoints
175 : INTEGER, DIMENSION(2, 3) :: bo
176 : REAL(kind=dp) :: a1, a2, epsilon_drho, epsilon_rho, gam, &
177 : sx
178 : REAL(kind=dp), CONTIGUOUS, DIMENSION(:, :, :), &
179 476 : POINTER :: e_0
180 5712 : TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_rho, ndrho, rho
181 : TYPE(xc_derivative_type), POINTER :: deriv
182 :
183 476 : NULLIFY (e_0)
184 1428 : DO ispin = 1, 2
185 952 : NULLIFY (e_rho(ispin)%array)
186 952 : NULLIFY (e_ndrho(ispin)%array)
187 952 : NULLIFY (rho(ispin)%array)
188 1428 : NULLIFY (ndrho(ispin)%array)
189 : END DO
190 :
191 476 : CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
192 476 : CALL section_vals_val_get(optx_params, "a1", r_val=a1)
193 476 : CALL section_vals_val_get(optx_params, "a2", r_val=a2)
194 476 : CALL section_vals_val_get(optx_params, "gamma", r_val=gam)
195 :
196 : CALL xc_rho_set_get(rho_set, rhoa=rho(1)%array, rhob=rho(2)%array, &
197 : norm_drhoa=ndrho(1)%array, &
198 : norm_drhob=ndrho(2)%array, rho_cutoff=epsilon_rho, &
199 476 : drho_cutoff=epsilon_drho, local_bounds=bo)
200 476 : npoints = (bo(2, 1) - bo(1, 1) + 1)*(bo(2, 2) - bo(1, 2) + 1)*(bo(2, 3) - bo(1, 3) + 1)
201 :
202 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], &
203 476 : allocate_deriv=.TRUE.)
204 476 : CALL xc_derivative_get(deriv, deriv_data=e_0)
205 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
206 476 : allocate_deriv=.TRUE.)
207 476 : CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
208 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
209 476 : allocate_deriv=.TRUE.)
210 476 : CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
211 :
212 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
213 476 : allocate_deriv=.TRUE.)
214 476 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
215 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
216 476 : allocate_deriv=.TRUE.)
217 476 : CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
218 :
219 476 : IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
220 0 : CPABORT("derivatives bigger than 1 not implemented")
221 : END IF
222 1428 : DO ispin = 1, 2
223 : CALL optx_lsd_calc(rho=rho(ispin)%array, norm_drho=ndrho(ispin)%array, &
224 : e_0=e_0, e_rho=e_rho(ispin)%array, e_ndrho=e_ndrho(ispin)%array, &
225 : npoints=npoints, epsilon_rho=epsilon_rho, &
226 : epsilon_drho=epsilon_drho, sx=sx, &
227 1428 : a1=a1, a2=a2, gam=gam)
228 : END DO
229 476 : END SUBROUTINE optx_lsd_eval
230 :
231 : ! **************************************************************************************************
232 : !> \brief optx exchange functional
233 : !> \param rho the full density
234 : !> \param norm_drho the norm of the gradient of the full density
235 : !> \param e_0 the value of the functional in that point
236 : !> \param e_rho the derivative of the functional wrt. rho
237 : !> \param e_ndrho the derivative of the functional wrt. norm_drho
238 : !> \param epsilon_rho the cutoff on rho
239 : !> \param epsilon_drho ...
240 : !> \param npoints ...
241 : !> \param sx scaling-parameter for exchange
242 : !> \param a1 a1 coefficient of the OPTX functional
243 : !> \param a2 a2 coefficient of the OPTX functional
244 : !> \param gam gamma coefficient of the OPTX functional
245 : !> \par History
246 : !> 01.2007 added scaling [Manuel Guidon]
247 : !> \author Joost VandeVondele
248 : ! **************************************************************************************************
249 564 : SUBROUTINE optx_lda_calc(rho, norm_drho, e_0, e_rho, e_ndrho, &
250 : epsilon_rho, epsilon_drho, npoints, sx, &
251 : a1, a2, gam)
252 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho
253 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho
254 : REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho
255 : INTEGER, INTENT(in) :: npoints
256 : REAL(kind=dp), INTENT(in) :: sx, a1, a2, gam
257 :
258 : REAL(KIND=dp), PARAMETER :: cx = 0.930525736349100_dp, &
259 : o43 = 4.0_dp/3.0_dp
260 :
261 : INTEGER :: ii
262 : REAL(KIND=dp) :: denom, ex, gamxsxs, myndrho, myrho, &
263 : rho43, tmp, xs
264 :
265 : !$OMP PARALLEL DO DEFAULT (NONE) &
266 : !$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
267 : !$OMP SHARED(epsilon_rho, epsilon_drho, sx, npoints) &
268 : !$OMP SHARED(a1, a2, gam) &
269 : !$OMP PRIVATE(ii, myrho, myndrho, rho43, xs, gamxsxs) &
270 564 : !$OMP PRIVATE(denom, ex, tmp)
271 :
272 : DO ii = 1, npoints
273 : ! we get the full density and need spin parts -> 0.5
274 : myrho = 0.5_dp*rho(ii)
275 : myndrho = 0.5_dp*MAX(norm_drho(ii), epsilon_drho)
276 : IF (myrho > 0.5_dp*epsilon_rho) THEN
277 : rho43 = myrho**o43
278 : xs = (myndrho/rho43)
279 : gamxsxs = gam*xs*xs
280 : denom = 1.0_dp/(1.0_dp + gamxsxs)
281 : ex = rho43*(a1*cx + a2*(gamxsxs*denom)**2)
282 : ! 2.0 for both spins
283 : e_0(ii) = e_0(ii) - (2.0_dp*ex)*sx
284 : tmp = rho43*2.0_dp*a2*gamxsxs*denom**2*(1.0_dp - gamxsxs*denom)
285 : ! derive e_0 wrt to rho (full) and ndrho (also full)
286 : e_rho(ii) = e_rho(ii) - ((o43*ex + tmp*gamxsxs*(-2.0_dp*o43))/myrho)*sx
287 : e_ndrho(ii) = e_ndrho(ii) - ((tmp*gam*2.0_dp*myndrho/rho43**2))*sx
288 : END IF
289 : END DO
290 :
291 : !$OMP END PARALLEL DO
292 :
293 564 : END SUBROUTINE optx_lda_calc
294 :
295 : ! **************************************************************************************************
296 : !> \brief optx exchange functional
297 : !> \param rho the *spin* density
298 : !> \param norm_drho the norm of the gradient of the *spin* density
299 : !> \param e_0 the value of the functional in that point
300 : !> \param e_rho the derivative of the functional wrt. rho
301 : !> \param e_ndrho the derivative of the functional wrt. norm_drho
302 : !> \param epsilon_rho the cutoff on rho
303 : !> \param epsilon_drho ...
304 : !> \param npoints ...
305 : !> \param sx scaling parameter for exchange
306 : !> \param a1 a1 coefficient of the OPTX functional
307 : !> \param a2 a2 coefficient of the OPTX functional
308 : !> \param gam gamma coefficient of the OPTX functional
309 : !> \par History
310 : !> 01.2007 added scaling [Manuel Guidon]
311 : !> \author Joost VandeVondele
312 : ! **************************************************************************************************
313 952 : SUBROUTINE optx_lsd_calc(rho, norm_drho, e_0, e_rho, e_ndrho, &
314 : epsilon_rho, epsilon_drho, npoints, sx, &
315 : a1, a2, gam)
316 : REAL(KIND=dp), DIMENSION(*), INTENT(IN) :: rho, norm_drho
317 : REAL(KIND=dp), DIMENSION(*), INTENT(INOUT) :: e_0, e_rho, e_ndrho
318 : REAL(kind=dp), INTENT(in) :: epsilon_rho, epsilon_drho
319 : INTEGER, INTENT(in) :: npoints
320 : REAL(kind=dp), INTENT(in) :: sx, a1, a2, gam
321 :
322 : REAL(KIND=dp), PARAMETER :: cx = 0.930525736349100_dp, &
323 : o43 = 4.0_dp/3.0_dp
324 :
325 : INTEGER :: ii
326 : REAL(KIND=dp) :: denom, ex, gamxsxs, myndrho, myrho, &
327 : rho43, tmp, xs
328 :
329 : !$OMP PARALLEL DO DEFAULT(NONE) &
330 : !$OMP SHARED(rho, norm_drho, e_0, e_rho, e_ndrho) &
331 : !$OMP SHARED(epsilon_rho, epsilon_drho, npoints, sx) &
332 : !$OMP SHARED(a1, a2, gam) &
333 : !$OMP PRIVATE(ii, denom, ex, gamxsxs, myndrho, myrho) &
334 952 : !$OMP PRIVATE(rho43, tmp, xs)
335 :
336 : DO ii = 1, npoints
337 : ! we do have the spin density already
338 : myrho = rho(ii)
339 : myndrho = MAX(norm_drho(ii), epsilon_drho)
340 : IF (myrho > epsilon_rho) THEN
341 : rho43 = myrho**o43
342 : xs = (myndrho/rho43)
343 : gamxsxs = gam*xs*xs
344 : denom = 1.0_dp/(1.0_dp + gamxsxs)
345 : ex = rho43*(a1*cx + a2*(gamxsxs*denom)**2)
346 : ! for a single spin
347 : e_0(ii) = e_0(ii) - ex*sx
348 : tmp = rho43*2.0_dp*a2*gamxsxs*denom**2*(1.0_dp - gamxsxs*denom)
349 : ! derive e_0 wrt to rho and ndrho
350 : e_rho(ii) = e_rho(ii) - ((o43*ex + tmp*gamxsxs*(-2.0_dp*o43))/myrho)*sx
351 : e_ndrho(ii) = e_ndrho(ii) - ((tmp*gam*2.0_dp*myndrho/rho43**2))*sx
352 : END IF
353 : END DO
354 :
355 : !$OMP END PARALLEL DO
356 :
357 952 : END SUBROUTINE optx_lsd_calc
358 :
359 : END MODULE xc_optx
|
