(git:e7e05ae)
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 ! **************************************************************************************************
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  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 
359 END 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