(git:374b731)
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-2024 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 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 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 POINTER :: e_0, e_ndrho, e_rho, norm_drho, rho
121 TYPE(xc_derivative_type), POINTER :: deriv
122
123 NULLIFY (e_0, e_ndrho, e_rho, norm_drho, rho)
124
125 CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
126 CALL section_vals_val_get(optx_params, "a1", r_val=a1)
127 CALL section_vals_val_get(optx_params, "a2", r_val=a2)
128 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 drho_cutoff=epsilon_drho)
133 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 allocate_deriv=.true.)
137 CALL xc_derivative_get(deriv, deriv_data=e_0)
138 deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], &
139 allocate_deriv=.true.)
140 CALL xc_derivative_get(deriv, deriv_data=e_rho)
141 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], &
142 allocate_deriv=.true.)
143 CALL xc_derivative_get(deriv, deriv_data=e_ndrho)
144 IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
145 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 a1=a1, a2=a2, gam=gam)
153 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 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 POINTER :: e_0
180 TYPE(cp_3d_r_cp_type), DIMENSION(2) :: e_ndrho, e_rho, ndrho, rho
181 TYPE(xc_derivative_type), POINTER :: deriv
182
183 NULLIFY (e_0)
184 DO ispin = 1, 2
185 NULLIFY (e_rho(ispin)%array)
186 NULLIFY (e_ndrho(ispin)%array)
187 NULLIFY (rho(ispin)%array)
188 NULLIFY (ndrho(ispin)%array)
189 END DO
190
191 CALL section_vals_val_get(optx_params, "scale_x", r_val=sx)
192 CALL section_vals_val_get(optx_params, "a1", r_val=a1)
193 CALL section_vals_val_get(optx_params, "a2", r_val=a2)
194 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 drho_cutoff=epsilon_drho, local_bounds=bo)
200 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 allocate_deriv=.true.)
204 CALL xc_derivative_get(deriv, deriv_data=e_0)
205 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], &
206 allocate_deriv=.true.)
207 CALL xc_derivative_get(deriv, deriv_data=e_rho(1)%array)
208 deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], &
209 allocate_deriv=.true.)
210 CALL xc_derivative_get(deriv, deriv_data=e_rho(2)%array)
211
212 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], &
213 allocate_deriv=.true.)
214 CALL xc_derivative_get(deriv, deriv_data=e_ndrho(1)%array)
215 deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], &
216 allocate_deriv=.true.)
217 CALL xc_derivative_get(deriv, deriv_data=e_ndrho(2)%array)
218
219 IF (grad_deriv > 1 .OR. grad_deriv < -1) THEN
220 cpabort("derivatives bigger than 1 not implemented")
221 END IF
222 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 a1=a1, a2=a2, gam=gam)
228 END DO
229 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 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!$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 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 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!$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 END SUBROUTINE optx_lsd_calc
358
359END 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:1040
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:66
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:169
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:110
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:126
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:49
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