(git:6344619)
Loading...
Searching...
No Matches
xc_optx.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE 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
40CONTAINS
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 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 IF (PRESENT(reference)) THEN
57 reference = "OPTX, Handy NC and Cohen AJ, JCP 116, p. 5411 (2002) (LDA)"
58 END IF
59 IF (PRESENT(shortform)) THEN
60 shortform = "OPTX exchange (LDA)"
61 END IF
62 IF (PRESENT(needs)) THEN
63 needs%rho = .true.
64 needs%norm_drho = .true.
65 END IF
66 IF (PRESENT(max_deriv)) max_deriv = 1
67 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 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 IF (PRESENT(reference)) THEN
84 reference = "OPTX, Handy NC and Cohen AJ, JCP 116, p. 5411 (2002), (LSD) "
85 END IF
86 IF (PRESENT(shortform)) THEN
87 shortform = "OPTX exchange (LSD)"
88 END IF
89 IF (PRESENT(needs)) THEN
90 needs%rho_spin = .true.
91 needs%norm_drho_spin = .true.
92 END IF
93 IF (PRESENT(max_deriv)) max_deriv = 1
94 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 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 POINTER :: e_0, e_ndrho, e_rho, norm_drho, rho
119 TYPE(xc_derivative_type), POINTER :: deriv
120
121 NULLIFY (e_0, e_ndrho, e_rho, norm_drho, rho)
122
123 CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
124 CALL section_vals_val_get(optx_params, "a1", r_val=a1)
125 CALL section_vals_val_get(optx_params, "a2", r_val=a2)
126 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 drho_cutoff=epsilon_drho)
131 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 allocate_deriv=.true.)
135 CALL xc_derivative_get(deriv, deriv_data=e_0)
136 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
137 allocate_deriv=.true.)
138 CALL xc_derivative_get(deriv, deriv_data=e_rho)
139 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
140 allocate_deriv=.true.)
141 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
142 IF (grad_deriv > 1) THEN
143 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 a1=a1, a2=a2, gam=gam)
151 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 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 POINTER :: e_0
176 TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_rho, ndrho, rho
177 TYPE(xc_derivative_type), POINTER :: deriv
178
179 NULLIFY (e_0)
180 DO ispin = 1, 2
181 NULLIFY (e_rho(ispin)%array)
182 NULLIFY (e_ndrho(ispin)%array)
183 NULLIFY (rho(ispin)%array)
184 NULLIFY (ndrho(ispin)%array)
185 END DO
186
187 CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
188 CALL section_vals_val_get(optx_params, "a1", r_val=a1)
189 CALL section_vals_val_get(optx_params, "a2", r_val=a2)
190 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 drho_cutoff=epsilon_drho, local_bounds=bo)
196 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 allocate_deriv=.true.)
200 CALL xc_derivative_get(deriv, deriv_data=e_0)
201 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
202 allocate_deriv=.true.)
203 CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
204 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
205 allocate_deriv=.true.)
206 CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
207
208 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
209 allocate_deriv=.true.)
210 CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
211 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
212 allocate_deriv=.true.)
213 CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
214
215 IF (grad_deriv > 1) THEN
216 cpabort("derivatives bigger than 1 not implemented")
217 END IF
218 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 a1=a1, a2=a2, gam=gam)
224 END DO
225 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 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!$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 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 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!$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 END SUBROUTINE optx_lsd_calc
354
355END MODULE xc_optx
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
xc_derivative_desc
Module with functions to handle derivative descriptors. derivative description are strings have the f...
Definition xc_derivative_desc.F:26
xc_derivative_desc::deriv_norm_drho
integer, parameter, public deriv_norm_drho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhoa
integer, parameter, public deriv_norm_drhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhob
integer, parameter, public deriv_rhob
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rhoa
integer, parameter, public deriv_rhoa
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_rho
integer, parameter, public deriv_rho
Definition xc_derivative_desc.F:35
xc_derivative_desc::deriv_norm_drhob
integer, parameter, public deriv_norm_drhob
Definition xc_derivative_desc.F:35
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_get_derivative
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
Definition xc_derivative_set_types.F:65
xc_derivative_types
Provides types for the management of the xc-functionals and their derivatives.
Definition xc_derivative_types.F:12
xc_derivative_types::xc_derivative_get
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
Definition xc_derivative_types.F:103
xc_optx
calculate optx
Definition xc_optx.F:14
xc_optx::optx_lda_info
subroutine, public optx_lda_info(reference, shortform, needs, max_deriv)
info about the optx functional
Definition xc_optx.F:52
xc_optx::optx_lsd_eval
subroutine, public optx_lsd_eval(rho_set, deriv_set, grad_deriv, optx_params)
evaluates the optx functional for lsd
Definition xc_optx.F:165
xc_optx::optx_lda_eval
subroutine, public optx_lda_eval(rho_set, deriv_set, grad_deriv, optx_params)
evaluates the optx functional for lda
Definition xc_optx.F:108
xc_optx::optx_lsd_info
subroutine, public optx_lsd_info(reference, shortform, needs, max_deriv)
info about the optx functional (LSD)
Definition xc_optx.F:79
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition xc_rho_set_types.F:281
cp_array_utils::cp_3d_r_cp_type
represent a pointer to a contiguous 3d array
Definition cp_array_utils.F:149
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_derivative_types::xc_derivative_type
represent a derivative of a functional
Definition xc_derivative_types.F:32
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78