(git:b279b6b)
atom_xc.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 routines that build the integrals of the Vxc potential calculated
10 !> for the atomic code
11 ! **************************************************************************************************
12 MODULE atom_xc
13 
14  USE atom_types, ONLY: atom_basis_type,&
15  atom_type,&
16  gth_pseudo,&
17  opmat_type,&
18  sgp_pseudo,&
19  upf_pseudo
20  USE atom_utils, ONLY: atom_core_density,&
21  atom_density,&
22  coulomb_potential_numeric,&
23  integrate_grid,&
24  numpot_matrix
25  USE cp_files, ONLY: close_file,&
26  get_unit_number,&
27  open_file
28  USE input_constants, ONLY: xc_none
29  USE input_section_types, ONLY: section_vals_get_subs_vals,&
30  section_vals_type,&
31  section_vals_val_get
32  USE kinds, ONLY: dp
33  USE mathconstants, ONLY: fourpi,&
34  pi
35  USE xc_atom, ONLY: xc_rho_set_atom_update
36  USE xc_derivative_desc, ONLY: &
37  deriv_norm_drho, deriv_norm_drhoa, deriv_norm_drhob, deriv_rho, deriv_rhoa, deriv_rhob, &
38  deriv_tau, deriv_tau_a, deriv_tau_b
39  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
40  xc_dset_create,&
41  xc_dset_get_derivative,&
42  xc_dset_release,&
43  xc_dset_zero_all
44  USE xc_derivative_types, ONLY: xc_derivative_get,&
45  xc_derivative_type
46  USE xc_derivatives, ONLY: xc_functionals_eval,&
47  xc_functionals_get_needs
48  USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
49  USE xc_rho_set_types, ONLY: xc_rho_set_create,&
50  xc_rho_set_release,&
51  xc_rho_set_type
52 #include "./base/base_uses.f90"
53 
54  IMPLICIT NONE
55 
56  PRIVATE
57 
58  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_xc'
59 
60 ! ZMP added coulomb_potential_numeric
61 ! ZMP added public subroutines calculate_atom_zmp, calculate_atom_ext_vxc
62  PUBLIC :: calculate_atom_vxc_lda, calculate_atom_vxc_lsd, &
63  calculate_atom_zmp, calculate_atom_ext_vxc
64  PUBLIC :: atom_vxc_lda, atom_vxc_lsd, atom_dpot_lda
65 
66 CONTAINS
67 
68 ! **************************************************************************************************
69 !> \brief ZMP subroutine for the calculation of the effective constraint
70 !> potential in the atomic code.
71 !> \param ext_density external electron density
72 !> \param atom information about the atomic kind
73 !> \param lprint save exchange-correlation potential to the file 'linear_potentials.dat'
74 !> \param xcmat exchange-correlation potential matrix
75 !> \author D. Varsano [daniele.varsano@nano.cnr.it]
76 ! **************************************************************************************************
77  SUBROUTINE calculate_atom_zmp(ext_density, atom, lprint, xcmat)
78  REAL(kind=dp), DIMENSION(:), INTENT(INOUT) :: ext_density
79  TYPE(atom_type), INTENT(INOUT) :: atom
80  LOGICAL :: lprint
81  TYPE(opmat_type), OPTIONAL, POINTER :: xcmat
82 
83  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_atom_zmp'
84 
85  INTEGER :: extunit, handle, ir, nr, z
86  REAL(kind=dp) :: int1, int2
87  REAL(kind=dp), DIMENSION(:), POINTER :: deltarho, rho_dum, vxc, vxc1, vxc2
88  REAL(kind=dp), DIMENSION(:, :), POINTER :: rho
89 
90  CALL timeset(routinen, handle)
91 
92  nr = atom%basis%grid%nr
93  z = atom%z
94  ALLOCATE (rho(nr, 1), vxc(nr), vxc1(nr), vxc2(nr), rho_dum(nr), deltarho(nr))
95  CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
96  !Vxc1
97  int1 = integrate_grid(ext_density, atom%basis%grid)
98  int1 = fourpi*int1
99  CALL coulomb_potential_numeric(vxc1, rho(:, 1), atom%basis%grid)
100 
101  vxc1 = -vxc1/z
102 
103  !Vxc2
104  rho_dum = rho(:, 1)*int1/z
105  deltarho = rho_dum - ext_density
106  int2 = integrate_grid(deltarho, atom%basis%grid)
107  CALL coulomb_potential_numeric(vxc2, deltarho, atom%basis%grid)
108 
109  vxc2 = vxc2*atom%lambda
110 
111  !Vxc
112  vxc = vxc1 + vxc2
113 
114  atom%energy%exc = fourpi*integrate_grid(vxc, rho(:, 1), atom%basis%grid)
115  atom%rho_diff_integral = fourpi*int2
116 
117  IF (lprint) THEN
118  extunit = get_unit_number()
119  CALL open_file(file_name="linear_potentials.dat", file_status="UNKNOWN", &
120  file_form="FORMATTED", file_action="WRITE", &
121  unit_number=extunit)
122 
123  WRITE (extunit, fmt='("#",T10,"R[bohr]",T36,"Rho[au]",T61,"DRho[au]",T86,"V_XC[au]",T111,"V_XC1[au]",T136,"V_XC2[au]")')
124  DO ir = 1, nr
125  WRITE (extunit, fmt='(T1,E24.15,T26,E24.15,T51,E24.15,T76,E24.15,T101,E24.15,T126,E24.15)') &
126  atom%basis%grid%rad(ir), rho(ir, 1), deltarho(ir), vxc(ir), vxc1(ir), vxc2(ir)
127  END DO
128  CALL close_file(unit_number=extunit)
129  END IF
130 
131  IF (PRESENT(xcmat)) CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)
132 
133  DEALLOCATE (rho, vxc, vxc1, vxc2, rho_dum, deltarho)
134 
135  CALL timestop(handle)
136 
137  END SUBROUTINE calculate_atom_zmp
138 
139 ! **************************************************************************************************
140 !> \brief ZMP subroutine for the scf external density from the external v_xc.
141 !> \param vxc exchange-correlation potential
142 !> \param atom information about the atomic kind
143 !> \param lprint save exchange-correlation potential to the file 'linear_potentials.dat'
144 !> \param xcmat exchange-correlation potential matrix
145 !> \author D. Varsano [daniele.varsano@nano.cnr.it]
146 ! **************************************************************************************************
147  SUBROUTINE calculate_atom_ext_vxc(vxc, atom, lprint, xcmat)
148  REAL(kind=dp), DIMENSION(:), INTENT(IN) :: vxc
149  TYPE(atom_type), INTENT(INOUT) :: atom
150  LOGICAL, INTENT(in) :: lprint
151  TYPE(opmat_type), INTENT(inout), OPTIONAL, POINTER :: xcmat
152 
153  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_atom_ext_vxc'
154 
155  INTEGER :: extunit, handle, ir, nr
156  REAL(kind=dp), DIMENSION(:, :), POINTER :: rho
157 
158  CALL timeset(routinen, handle)
159  nr = atom%basis%grid%nr
160 
161  ALLOCATE (rho(nr, 1))
162 
163  CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
164 
165  IF (lprint) THEN
166  extunit = get_unit_number()
167  CALL open_file(file_name="linear_potentials.dat", file_status="UNKNOWN", &
168  file_form="FORMATTED", file_action="WRITE", &
169  unit_number=extunit)
170 
171  WRITE (extunit, fmt='("#",T10,"R[bohr]",T36,"Rho[au]",T61,"V_XC[au]")')
172  DO ir = 1, nr
173  WRITE (extunit, fmt='(T1,E24.15,T26,E24.15,T51,E24.15)') &
174  atom%basis%grid%rad(ir), rho(ir, 1), vxc(ir)
175  END DO
176  CALL close_file(unit_number=extunit)
177  END IF
178 
179  atom%energy%exc = fourpi*integrate_grid(vxc, rho(:, 1), atom%basis%grid)
180  IF (PRESENT(xcmat)) CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)
181 
182  DEALLOCATE (rho)
183 
184  CALL timestop(handle)
185 
186  END SUBROUTINE calculate_atom_ext_vxc
187 
188 ! **************************************************************************************************
189 !> \brief Calculate atomic exchange-correlation potential in spin-restricted case.
190 !> \param xcmat exchange-correlation potential expressed in the atomic basis set
191 !> \param atom information about the atomic kind
192 !> \param xc_section XC input section
193 !> \par History
194 !> * 03.2010 extension of GTH pseudo-potential definition [Juerg Hutter]
195 !> * 02.2010 renamed to calculate_atom_vxc_lda() [Juerg Hutter]
196 !> * 08.2008 created [Juerg Hutter]
197 ! **************************************************************************************************
198  SUBROUTINE calculate_atom_vxc_lda(xcmat, atom, xc_section)
199  TYPE(opmat_type), POINTER :: xcmat
200  TYPE(atom_type), INTENT(INOUT) :: atom
201  TYPE(section_vals_type), POINTER :: xc_section
202 
203  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_atom_vxc_lda'
204 
205  INTEGER :: deriv_order, handle, i, l, myfun, n1, &
206  n2, n3, nr, nspins, unit_nr
207  INTEGER, DIMENSION(2, 3) :: bounds
208  LOGICAL :: lsd, nlcc
209  REAL(kind=dp) :: density_cut, gradient_cut, tau_cut
210  REAL(kind=dp), DIMENSION(:), POINTER :: exc, vxc
211  REAL(kind=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
212  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: taumat, xcpot
213  TYPE(section_vals_type), POINTER :: xc_fun_section
214  TYPE(xc_derivative_set_type) :: deriv_set
215  TYPE(xc_derivative_type), POINTER :: deriv
216  TYPE(xc_rho_cflags_type) :: needs
217  TYPE(xc_rho_set_type) :: rho_set
218 
219  CALL timeset(routinen, handle)
220 
221  nlcc = .false.
222  IF (atom%potential%ppot_type == gth_pseudo) THEN
223  nlcc = atom%potential%gth_pot%nlcc
224  ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
225  nlcc = atom%potential%upf_pot%core_correction
226  ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
227  nlcc = atom%potential%sgp_pot%has_nlcc
228  END IF
229 
230  IF (ASSOCIATED(xc_section)) THEN
231  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
232  CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
233 
234  IF (myfun == xc_none) THEN
235  atom%energy%exc = 0._dp
236  ELSE
237  CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
238  CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
239  CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
240 
241  lsd = .false.
242  nspins = 1
243  needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.false.)
244 
245  ! Prepare the structures needed to calculate and store the xc derivatives
246 
247  ! Array dimension: here anly one dimensional arrays are used,
248  ! i.e. only the first column of deriv_data is read.
249  ! The other to dimensions are set to size equal 1
250  nr = atom%basis%grid%nr
251  bounds(1:2, 1:3) = 1
252  bounds(2, 1) = nr
253 
254  ! create a place where to put the derivatives
255  CALL xc_dset_create(deriv_set, local_bounds=bounds)
256  ! create the place where to store the argument for the functionals
257  CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
258  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
259  ! allocate the required 3d arrays where to store rho and drho
260  CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
261 
262  NULLIFY (rho, drho, tau)
263  IF (needs%rho) THEN
264  ALLOCATE (rho(nr, 1))
265  CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
266  IF (nlcc) THEN
267  CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
268  END IF
269  END IF
270  IF (needs%norm_drho) THEN
271  ALLOCATE (drho(nr, 1))
272  CALL atom_density(drho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="DER")
273  IF (nlcc) THEN
274  CALL atom_core_density(drho(:, 1), atom%potential, typ="DER", rr=atom%basis%grid%rad)
275  END IF
276  END IF
277  IF (needs%tau) THEN
278  ALLOCATE (tau(nr, 1))
279  CALL atom_density(tau(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
280  typ="KIN", rr=atom%basis%grid%rad2)
281  END IF
282  IF (needs%laplace_rho) THEN
283  ALLOCATE (lap(nr, 1))
284  CALL atom_density(lap(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
285  typ="LAP", rr=atom%basis%grid%rad2)
286  IF (nlcc) THEN
287  CALL atom_core_density(lap(:, 1), atom%potential, typ="LAP", rr=atom%basis%grid%rad)
288  END IF
289  END IF
290 
291  CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
292 
293  CALL xc_dset_zero_all(deriv_set)
294 
295  deriv_order = 1
296  CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
297  deriv_order=deriv_order)
298 
299  ! Integration to get the matrix elements and energy
300  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.false.)
301  CALL xc_derivative_get(deriv, deriv_data=xcpot)
302  atom%energy%exc = fourpi*integrate_grid(xcpot(:, 1, 1), atom%basis%grid)
303 
304  ! dump grid density and xcpot (xc energy?)
305  IF (.false.) THEN
306  CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atom.dat")
307  DO i = 1, SIZE(xcpot(:, 1, 1))
308  WRITE (unit_nr, *) atom%basis%grid%rad(i), rho(i, 1), xcpot(i, 1, 1)
309  END DO
310  CALL close_file(unit_nr)
311  END IF
312 
313  IF (needs%rho) THEN
314  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], allocate_deriv=.false.)
315  CALL xc_derivative_get(deriv, deriv_data=xcpot)
316  CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 0)
317  IF (.false.) THEN
318  CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atompot.dat")
319  DO i = 1, SIZE(xcpot(:, 1, 1))
320  WRITE (unit_nr, *) atom%basis%grid%rad(i), rho(i, 1), xcpot(i, 1, 1)
321  END DO
322  CALL close_file(unit_nr)
323  END IF
324  DEALLOCATE (rho)
325  END IF
326  IF (needs%norm_drho) THEN
327  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], allocate_deriv=.false.)
328  CALL xc_derivative_get(deriv, deriv_data=xcpot)
329  CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 1)
330  IF (.false.) THEN
331  CALL open_file(unit_number=unit_nr, file_status="UNKNOWN", file_action="WRITE", file_name="atomdpot.dat")
332  DO i = 1, SIZE(xcpot(:, 1, 1))
333  WRITE (unit_nr, *) atom%basis%grid%rad(i), drho(i, 1), xcpot(i, 1, 1)
334  END DO
335  CALL close_file(unit_nr)
336  END IF
337  DEALLOCATE (drho)
338  END IF
339  IF (needs%tau) THEN
340  deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], allocate_deriv=.false.)
341  CALL xc_derivative_get(deriv, deriv_data=xcpot)
342  n1 = SIZE(xcmat%op, 1)
343  n2 = SIZE(xcmat%op, 2)
344  n3 = SIZE(xcmat%op, 3)
345  ALLOCATE (taumat(n1, n2, 0:n3 - 1))
346  taumat = 0._dp
347 
348  xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
349  CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 2)
350  xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
351  CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
352  DO l = 0, 3
353  xcmat%op(:, :, l) = xcmat%op(:, :, l) + real(l*(l + 1), dp)*taumat(:, :, l)
354  END DO
355 
356  DEALLOCATE (tau)
357  DEALLOCATE (taumat)
358  END IF
359 
360  ! assume lap is not really needed
361  IF (needs%laplace_rho) THEN
362  DEALLOCATE (lap)
363  END IF
364 
365  ! Release the xc structure used to store the xc derivatives
366  CALL xc_dset_release(deriv_set)
367  CALL xc_rho_set_release(rho_set)
368 
369  END IF !xc_none
370 
371  ELSE
372 
373  ! we don't have an xc_section, use a default setup
374  nr = atom%basis%grid%nr
375  ALLOCATE (rho(nr, 1), exc(nr), vxc(nr))
376 
377  CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
378  IF (nlcc) THEN
379  CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
380  END IF
381  CALL lda_pade(rho(:, 1), exc, vxc)
382 
383  atom%energy%exc = fourpi*integrate_grid(exc, atom%basis%grid)
384  CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)
385 
386  DEALLOCATE (rho, exc, vxc)
387 
388  END IF
389 
390  CALL timestop(handle)
391 
392  END SUBROUTINE calculate_atom_vxc_lda
393 
394 ! **************************************************************************************************
395 !> \brief Calculate atomic exchange-correlation potential in spin-polarized case.
396 !> \param xcmata spin-alpha exchange-correlation potential matrix
397 !> \param xcmatb spin-beta exchange-correlation potential matrix
398 !> \param atom information about the atomic kind
399 !> \param xc_section XC input section
400 !> \par History
401 !> * 02.2010 created [Juerg Hutter]
402 ! **************************************************************************************************
403  SUBROUTINE calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
404  TYPE(opmat_type), POINTER :: xcmata, xcmatb
405  TYPE(atom_type), INTENT(INOUT) :: atom
406  TYPE(section_vals_type), POINTER :: xc_section
407 
408  CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_atom_vxc_lsd'
409 
410  INTEGER :: deriv_order, handle, l, myfun, n1, n2, &
411  n3, nr, nspins
412  INTEGER, DIMENSION(2, 3) :: bounds
413  LOGICAL :: lsd, nlcc
414  REAL(kind=dp) :: density_cut, gradient_cut, tau_cut
415  REAL(kind=dp), DIMENSION(:), POINTER :: exc, vxca, vxcb, xfun
416  REAL(kind=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
417  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: taumat, xcpot
418  TYPE(section_vals_type), POINTER :: xc_fun_section
419  TYPE(xc_derivative_set_type) :: deriv_set
420  TYPE(xc_derivative_type), POINTER :: deriv
421  TYPE(xc_rho_cflags_type) :: needs
422  TYPE(xc_rho_set_type) :: rho_set
423 
424  CALL timeset(routinen, handle)
425 
426  nlcc = .false.
427  IF (atom%potential%ppot_type == gth_pseudo) THEN
428  nlcc = atom%potential%gth_pot%nlcc
429  ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
430  nlcc = atom%potential%upf_pot%core_correction
431  cpassert(.NOT. nlcc)
432  ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
433  nlcc = atom%potential%sgp_pot%has_nlcc
434  cpassert(.NOT. nlcc)
435  END IF
436 
437  IF (ASSOCIATED(xc_section)) THEN
438  IF (nlcc) THEN
439  nr = atom%basis%grid%nr
440  ALLOCATE (xfun(nr))
441  END IF
442 
443  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
444  CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
445 
446  IF (myfun == xc_none) THEN
447  atom%energy%exc = 0._dp
448  ELSE
449  CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
450  CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
451  CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
452 
453  lsd = .true.
454  nspins = 2
455  needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.false.)
456 
457  ! Prepare the structures needed to calculate and store the xc derivatives
458 
459  ! Array dimension: here anly one dimensional arrays are used,
460  ! i.e. only the first column of deriv_data is read.
461  ! The other to dimensions are set to size equal 1
462  nr = atom%basis%grid%nr
463  bounds(1:2, 1:3) = 1
464  bounds(2, 1) = nr
465 
466  ! create a place where to put the derivatives
467  CALL xc_dset_create(deriv_set, local_bounds=bounds)
468  ! create the place where to store the argument for the functionals
469  CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
470  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
471  ! allocate the required 3d arrays where to store rho and drho
472  CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
473 
474  NULLIFY (rho, drho, tau)
475  IF (needs%rho_spin) THEN
476  ALLOCATE (rho(nr, 2))
477  CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
478  CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
479  IF (nlcc) THEN
480  xfun = 0.0_dp
481  CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
482  rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
483  rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
484  END IF
485  END IF
486  IF (needs%norm_drho_spin) THEN
487  ALLOCATE (drho(nr, 2))
488  CALL atom_density(drho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="DER")
489  CALL atom_density(drho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="DER")
490  IF (nlcc) THEN
491  xfun = 0.0_dp
492  CALL atom_core_density(xfun(:), atom%potential, typ="DER", rr=atom%basis%grid%rad)
493  drho(:, 1) = drho(:, 1) + 0.5_dp*xfun(:)
494  drho(:, 2) = drho(:, 2) + 0.5_dp*xfun(:)
495  END IF
496  END IF
497  IF (needs%tau_spin) THEN
498  ALLOCATE (tau(nr, 2))
499  CALL atom_density(tau(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, &
500  typ="KIN", rr=atom%basis%grid%rad2)
501  CALL atom_density(tau(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, &
502  typ="KIN", rr=atom%basis%grid%rad2)
503  END IF
504  IF (needs%laplace_rho_spin) THEN
505  ALLOCATE (lap(nr, 2))
506  CALL atom_density(lap(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, &
507  typ="LAP", rr=atom%basis%grid%rad2)
508  CALL atom_density(lap(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, &
509  typ="LAP", rr=atom%basis%grid%rad2)
510  IF (nlcc) THEN
511  xfun = 0.0_dp
512  CALL atom_core_density(xfun(:), atom%potential, typ="LAP", rr=atom%basis%grid%rad)
513  lap(:, 1) = lap(:, 1) + 0.5_dp*xfun(:)
514  lap(:, 2) = lap(:, 2) + 0.5_dp*xfun(:)
515  END IF
516  END IF
517 
518  CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
519 
520  CALL xc_dset_zero_all(deriv_set)
521 
522  deriv_order = 1
523  CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
524  deriv_order=deriv_order)
525 
526  ! Integration to get the matrix elements and energy
527  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.false.)
528  CALL xc_derivative_get(deriv, deriv_data=xcpot)
529  atom%energy%exc = fourpi*integrate_grid(xcpot(:, 1, 1), atom%basis%grid)
530 
531  IF (needs%rho_spin) THEN
532  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], allocate_deriv=.false.)
533  CALL xc_derivative_get(deriv, deriv_data=xcpot)
534  CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 0)
535  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], allocate_deriv=.false.)
536  CALL xc_derivative_get(deriv, deriv_data=xcpot)
537  CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 0)
538  DEALLOCATE (rho)
539  END IF
540  IF (needs%norm_drho_spin) THEN
541  ! drhoa
542  NULLIFY (deriv)
543  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], allocate_deriv=.false.)
544  CALL xc_derivative_get(deriv, deriv_data=xcpot)
545  CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 1)
546  ! drhob
547  NULLIFY (deriv)
548  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], allocate_deriv=.false.)
549  CALL xc_derivative_get(deriv, deriv_data=xcpot)
550  CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 1)
551  ! Cross Terms
552  NULLIFY (deriv)
553  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
554  IF (ASSOCIATED(deriv)) THEN
555  CALL xc_derivative_get(deriv, deriv_data=xcpot)
556  CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 1)
557  CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 1)
558  END IF
559  DEALLOCATE (drho)
560  END IF
561  IF (needs%tau_spin) THEN
562  n1 = SIZE(xcmata%op, 1)
563  n2 = SIZE(xcmata%op, 2)
564  n3 = SIZE(xcmata%op, 3)
565  ALLOCATE (taumat(n1, n2, 0:n3 - 1))
566 
567  deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], allocate_deriv=.false.)
568  CALL xc_derivative_get(deriv, deriv_data=xcpot)
569  taumat = 0._dp
570  xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
571  CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 2)
572  xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
573  CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
574  DO l = 0, 3
575  xcmata%op(:, :, l) = xcmata%op(:, :, l) + real(l*(l + 1), dp)*taumat(:, :, l)
576  END DO
577 
578  deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], allocate_deriv=.false.)
579  CALL xc_derivative_get(deriv, deriv_data=xcpot)
580  taumat = 0._dp
581  xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
582  CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 2)
583  xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
584  CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
585  DO l = 0, 3
586  xcmatb%op(:, :, l) = xcmatb%op(:, :, l) + real(l*(l + 1), dp)*taumat(:, :, l)
587  END DO
588 
589  DEALLOCATE (tau)
590  DEALLOCATE (taumat)
591  END IF
592 
593  ! assume lap is not really needed
594  IF (needs%laplace_rho_spin) THEN
595  DEALLOCATE (lap)
596  END IF
597 
598  ! Release the xc structure used to store the xc derivatives
599  CALL xc_dset_release(deriv_set)
600  CALL xc_rho_set_release(rho_set)
601 
602  END IF !xc_none
603 
604  IF (nlcc) DEALLOCATE (xfun)
605 
606  ELSE
607 
608  ! we don't have an xc_section, use a default setup
609  nr = atom%basis%grid%nr
610  ALLOCATE (rho(nr, 2), exc(nr), vxca(nr), vxcb(nr))
611  IF (nlcc) ALLOCATE (xfun(nr))
612 
613  CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
614  CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
615  IF (nlcc) THEN
616  xfun(:) = 0.0_dp
617  CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
618  rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
619  rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
620  END IF
621  CALL lsd_pade(rho(:, 1), rho(:, 2), exc, vxca, vxcb)
622 
623  atom%energy%exc = fourpi*integrate_grid(exc, atom%basis%grid)
624  CALL numpot_matrix(xcmata%op, vxca, atom%basis, 0)
625  CALL numpot_matrix(xcmatb%op, vxcb, atom%basis, 0)
626 
627  DEALLOCATE (rho, exc, vxca, vxcb)
628  IF (nlcc) DEALLOCATE (xfun)
629 
630  END IF
631 
632  CALL timestop(handle)
633 
634  END SUBROUTINE calculate_atom_vxc_lsd
635 
636 ! **************************************************************************************************
637 !> \brief Alternative subroutine that calculates atomic exchange-correlation potential
638 !> in spin-restricted case.
639 !> \param basis atomic basis set
640 !> \param pmat density matrix
641 !> \param maxl maximum angular momentum
642 !> \param xc_section XC input section
643 !> \param fexc exchange-correlation energy
644 !> \param xcmat exchange-correlation potential matrix
645 !> \par History
646 !> * 07.2016 ADMM analysis [Juerg Hutter]
647 ! **************************************************************************************************
648  SUBROUTINE atom_vxc_lda(basis, pmat, maxl, xc_section, fexc, xcmat)
649  TYPE(atom_basis_type), POINTER :: basis
650  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: pmat
651  INTEGER, INTENT(IN) :: maxl
652  TYPE(section_vals_type), POINTER :: xc_section
653  REAL(kind=dp), INTENT(OUT) :: fexc
654  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: xcmat
655 
656  CHARACTER(LEN=*), PARAMETER :: routinen = 'atom_vxc_lda'
657 
658  INTEGER :: deriv_order, handle, l, myfun, n1, n2, &
659  n3, nr, nspins
660  INTEGER, DIMENSION(2, 3) :: bounds
661  LOGICAL :: lsd
662  REAL(kind=dp) :: density_cut, gradient_cut, tau_cut
663  REAL(kind=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
664  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: taumat, xcpot
665  TYPE(section_vals_type), POINTER :: xc_fun_section
666  TYPE(xc_derivative_set_type) :: deriv_set
667  TYPE(xc_derivative_type), POINTER :: deriv
668  TYPE(xc_rho_cflags_type) :: needs
669  TYPE(xc_rho_set_type) :: rho_set
670 
671  CALL timeset(routinen, handle)
672 
673  cpassert(ASSOCIATED(xc_section))
674 
675  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
676  CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
677 
678  IF (myfun == xc_none) THEN
679  fexc = 0._dp
680  ELSE
681  CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
682  CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
683  CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
684 
685  lsd = .false.
686  nspins = 1
687  needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.false.)
688 
689  ! Prepare the structures needed to calculate and store the xc derivatives
690 
691  ! Array dimension: here anly one dimensional arrays are used,
692  ! i.e. only the first column of deriv_data is read.
693  ! The other to dimensions are set to size equal 1
694  nr = basis%grid%nr
695  bounds(1:2, 1:3) = 1
696  bounds(2, 1) = nr
697 
698  ! create a place where to put the derivatives
699  CALL xc_dset_create(deriv_set, local_bounds=bounds)
700  ! create the place where to store the argument for the functionals
701  CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
702  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
703  ! allocate the required 3d arrays where to store rho and drho
704  CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
705 
706  NULLIFY (rho, drho, tau)
707  IF (needs%rho) THEN
708  ALLOCATE (rho(nr, 1))
709  CALL atom_density(rho(:, 1), pmat, basis, maxl, typ="RHO")
710  END IF
711  IF (needs%norm_drho) THEN
712  ALLOCATE (drho(nr, 1))
713  CALL atom_density(drho(:, 1), pmat, basis, maxl, typ="DER")
714  END IF
715  IF (needs%tau) THEN
716  ALLOCATE (tau(nr, 1))
717  CALL atom_density(tau(:, 1), pmat, basis, maxl, typ="KIN", rr=basis%grid%rad2)
718  END IF
719  IF (needs%laplace_rho) THEN
720  ALLOCATE (lap(nr, 1))
721  CALL atom_density(lap(:, 1), pmat, basis, maxl, typ="LAP", rr=basis%grid%rad2)
722  END IF
723 
724  CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
725 
726  CALL xc_dset_zero_all(deriv_set)
727 
728  deriv_order = 1
729  CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
730  deriv_order=deriv_order)
731 
732  ! Integration to get the matrix elements and energy
733  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.false.)
734  CALL xc_derivative_get(deriv, deriv_data=xcpot)
735  fexc = fourpi*integrate_grid(xcpot(:, 1, 1), basis%grid)
736 
737  IF (needs%rho) THEN
738  deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], allocate_deriv=.false.)
739  CALL xc_derivative_get(deriv, deriv_data=xcpot)
740  CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 0)
741  DEALLOCATE (rho)
742  END IF
743  IF (needs%norm_drho) THEN
744  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], allocate_deriv=.false.)
745  CALL xc_derivative_get(deriv, deriv_data=xcpot)
746  CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 1)
747  DEALLOCATE (drho)
748  END IF
749  IF (needs%tau) THEN
750  deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], allocate_deriv=.false.)
751  CALL xc_derivative_get(deriv, deriv_data=xcpot)
752  n1 = SIZE(xcmat, 1)
753  n2 = SIZE(xcmat, 2)
754  n3 = SIZE(xcmat, 3)
755  ALLOCATE (taumat(n1, n2, 0:n3 - 1))
756  taumat = 0._dp
757 
758  xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
759  CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 2)
760  xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
761  CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
762  DO l = 0, 3
763  xcmat(:, :, l) = xcmat(:, :, l) + real(l*(l + 1), dp)*taumat(:, :, l)
764  END DO
765 
766  DEALLOCATE (tau)
767  DEALLOCATE (taumat)
768  END IF
769 
770  ! we assume lap is not really needed
771  IF (needs%laplace_rho) THEN
772  DEALLOCATE (lap)
773  END IF
774 
775  ! Release the xc structure used to store the xc derivatives
776  CALL xc_dset_release(deriv_set)
777  CALL xc_rho_set_release(rho_set)
778 
779  END IF !xc_none
780 
781  CALL timestop(handle)
782 
783  END SUBROUTINE atom_vxc_lda
784 
785 ! **************************************************************************************************
786 !> \brief Alternative subroutine that calculates atomic exchange-correlation potential
787 !> in spin-polarized case.
788 !> \param basis atomic basis set
789 !> \param pmata spin-alpha density matrix
790 !> \param pmatb spin-beta density matrix
791 !> \param maxl maximum angular momentum
792 !> \param xc_section XC input section
793 !> \param fexc exchange-correlation energy
794 !> \param xcmata spin-alpha exchange-correlation potential matrix
795 !> \param xcmatb spin-beta exchange-correlation potential matrix
796 !> \par History
797 !> * 07.2016 ADMM analysis [Juerg Hutter]
798 ! **************************************************************************************************
799  SUBROUTINE atom_vxc_lsd(basis, pmata, pmatb, maxl, xc_section, fexc, xcmata, xcmatb)
800  TYPE(atom_basis_type), POINTER :: basis
801  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: pmata, pmatb
802  INTEGER, INTENT(IN) :: maxl
803  TYPE(section_vals_type), POINTER :: xc_section
804  REAL(kind=dp), INTENT(OUT) :: fexc
805  REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: xcmata, xcmatb
806 
807  CHARACTER(LEN=*), PARAMETER :: routinen = 'atom_vxc_lsd'
808 
809  INTEGER :: deriv_order, handle, l, myfun, n1, n2, &
810  n3, nr, nspins
811  INTEGER, DIMENSION(2, 3) :: bounds
812  LOGICAL :: lsd
813  REAL(kind=dp) :: density_cut, gradient_cut, tau_cut
814  REAL(kind=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
815  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: taumat, xcpot
816  TYPE(section_vals_type), POINTER :: xc_fun_section
817  TYPE(xc_derivative_set_type) :: deriv_set
818  TYPE(xc_derivative_type), POINTER :: deriv
819  TYPE(xc_rho_cflags_type) :: needs
820  TYPE(xc_rho_set_type) :: rho_set
821 
822  CALL timeset(routinen, handle)
823 
824  cpassert(ASSOCIATED(xc_section))
825 
826  xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
827  CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
828 
829  IF (myfun == xc_none) THEN
830  fexc = 0._dp
831  ELSE
832  CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
833  CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
834  CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
835 
836  lsd = .true.
837  nspins = 2
838  needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.false.)
839 
840  ! Prepare the structures needed to calculate and store the xc derivatives
841 
842  ! Array dimension: here anly one dimensional arrays are used,
843  ! i.e. only the first column of deriv_data is read.
844  ! The other to dimensions are set to size equal 1
845  nr = basis%grid%nr
846  bounds(1:2, 1:3) = 1
847  bounds(2, 1) = nr
848 
849  ! create a place where to put the derivatives
850  CALL xc_dset_create(deriv_set, local_bounds=bounds)
851  ! create the place where to store the argument for the functionals
852  CALL xc_rho_set_create(rho_set, bounds, rho_cutoff=density_cut, &
853  drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
854  ! allocate the required 3d arrays where to store rho and drho
855  CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
856 
857  NULLIFY (rho, drho, tau)
858  IF (needs%rho_spin) THEN
859  ALLOCATE (rho(nr, 2))
860  CALL atom_density(rho(:, 1), pmata, basis, maxl, typ="RHO")
861  CALL atom_density(rho(:, 2), pmatb, basis, maxl, typ="RHO")
862  END IF
863  IF (needs%norm_drho_spin) THEN
864  ALLOCATE (drho(nr, 2))
865  CALL atom_density(drho(:, 1), pmata, basis, maxl, typ="DER")
866  CALL atom_density(drho(:, 2), pmatb, basis, maxl, typ="DER")
867  END IF
868  IF (needs%tau_spin) THEN
869  ALLOCATE (tau(nr, 2))
870  CALL atom_density(tau(:, 1), pmata, basis, maxl, typ="KIN", rr=basis%grid%rad2)
871  CALL atom_density(tau(:, 2), pmatb, basis, maxl, typ="KIN", rr=basis%grid%rad2)
872  END IF
873  IF (needs%laplace_rho_spin) THEN
874  ALLOCATE (lap(nr, 2))
875  CALL atom_density(lap(:, 1), pmata, basis, maxl, typ="LAP", rr=basis%grid%rad2)
876  CALL atom_density(lap(:, 2), pmatb, basis, maxl, typ="LAP", rr=basis%grid%rad2)
877  END IF
878 
879  CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
880 
881  CALL xc_dset_zero_all(deriv_set)
882 
883  deriv_order = 1
884  CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
885  deriv_order=deriv_order)
886 
887  ! Integration to get the matrix elements and energy
888  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.false.)
889  CALL xc_derivative_get(deriv, deriv_data=xcpot)
890  fexc = fourpi*integrate_grid(xcpot(:, 1, 1), basis%grid)
891 
892  IF (needs%rho_spin) THEN
893  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], allocate_deriv=.false.)
894  CALL xc_derivative_get(deriv, deriv_data=xcpot)
895  CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 0)
896  deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], allocate_deriv=.false.)
897  CALL xc_derivative_get(deriv, deriv_data=xcpot)
898  CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 0)
899  DEALLOCATE (rho)
900  END IF
901  IF (needs%norm_drho_spin) THEN
902  ! drhoa
903  NULLIFY (deriv)
904  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], allocate_deriv=.false.)
905  CALL xc_derivative_get(deriv, deriv_data=xcpot)
906  CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 1)
907  ! drhob
908  NULLIFY (deriv)
909  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], allocate_deriv=.false.)
910  CALL xc_derivative_get(deriv, deriv_data=xcpot)
911  CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 1)
912  ! Cross Terms
913  NULLIFY (deriv)
914  deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
915  IF (ASSOCIATED(deriv)) THEN
916  CALL xc_derivative_get(deriv, deriv_data=xcpot)
917  CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 1)
918  CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 1)
919  END IF
920  DEALLOCATE (drho)
921  END IF
922  IF (needs%tau_spin) THEN
923  n1 = SIZE(xcmata, 1)
924  n2 = SIZE(xcmata, 2)
925  n3 = SIZE(xcmata, 3)
926  ALLOCATE (taumat(n1, n2, 0:n3 - 1))
927 
928  deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], allocate_deriv=.false.)
929  CALL xc_derivative_get(deriv, deriv_data=xcpot)
930  taumat = 0._dp
931  xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
932  CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 2)
933  xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
934  CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
935  DO l = 0, 3
936  xcmata(:, :, l) = xcmata(:, :, l) + real(l*(l + 1), dp)*taumat(:, :, l)
937  END DO
938 
939  deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], allocate_deriv=.false.)
940  CALL xc_derivative_get(deriv, deriv_data=xcpot)
941  taumat = 0._dp
942  xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
943  CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 2)
944  xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
945  CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
946  DO l = 0, 3
947  xcmatb(:, :, l) = xcmatb(:, :, l) + real(l*(l + 1), dp)*taumat(:, :, l)
948  END DO
949 
950  DEALLOCATE (tau)
951  DEALLOCATE (taumat)
952  END IF
953 
954  ! Release the xc structure used to store the xc derivatives
955  CALL xc_dset_release(deriv_set)
956  CALL xc_rho_set_release(rho_set)
957 
958  END IF !xc_none
959 
960  CALL timestop(handle)
961 
962  END SUBROUTINE atom_vxc_lsd
963 
964 ! **************************************************************************************************
965 !> \brief Estimate the ADMM exchange energy correction term using a model exchange functional.
966 !> \param basis0 reference atomic basis set
967 !> \param pmat0 reference density matrix
968 !> \param basis1 ADMM basis set
969 !> \param pmat1 ADMM density matrix
970 !> \param maxl maxumum angular momentum
971 !> \param functional name of the model exchange functional:
972 !> "LINX" -- linear extrapolation of the Slater exchange potential [?]
973 !> \param dfexc exchange-correlation energy difference
974 !> \param linxpar LINX parameters
975 !> \par History
976 !> * 07.2016 ADMM analysis [Juerg Hutter]
977 ! **************************************************************************************************
978  SUBROUTINE atom_dpot_lda(basis0, pmat0, basis1, pmat1, maxl, functional, dfexc, linxpar)
979  TYPE(atom_basis_type), POINTER :: basis0
980  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: pmat0
981  TYPE(atom_basis_type), POINTER :: basis1
982  REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: pmat1
983  INTEGER, INTENT(IN) :: maxl
984  CHARACTER(LEN=*), INTENT(IN) :: functional
985  REAL(kind=dp), INTENT(OUT) :: dfexc
986  REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: linxpar
987 
988  CHARACTER(LEN=*), PARAMETER :: routinen = 'atom_dpot_lda'
989  REAL(kind=dp), PARAMETER :: alx = 0.20_dp, blx = -0.06_dp, &
990  fs = -0.75_dp*(3._dp/pi)**(1._dp/3._dp)
991 
992  INTEGER :: handle, ir, nr
993  REAL(kind=dp) :: a, b, fx
994  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: delta, drho0, drho1, pot0, pot1, rho0, &
995  rho1
996 
997  CALL timeset(routinen, handle)
998 
999  nr = basis0%grid%nr
1000 
1001  ALLOCATE (rho0(nr), drho0(nr))
1002  CALL atom_density(rho0, pmat0, basis0, maxl, typ="RHO")
1003  CALL atom_density(drho0, pmat0, basis0, maxl, typ="DER")
1004  !
1005  ALLOCATE (rho1(nr), drho1(nr))
1006  CALL atom_density(rho1, pmat1, basis1, maxl, typ="RHO")
1007  ! CONSIDER TO REMOVE [?]: drho1 is calculated but it is not used.
1008  CALL atom_density(drho1, pmat1, basis1, maxl, typ="DER")
1009  !
1010  ALLOCATE (delta(nr))
1011  fx = 4.0_dp/3.0_dp
1012  delta(1:nr) = fs*(rho0(1:nr)**fx - rho1(1:nr)**fx)
1013 
1014  SELECT CASE (trim(functional))
1015  CASE ("LINX")
1016  IF (PRESENT(linxpar)) THEN
1017  a = linxpar(1)
1018  b = linxpar(2)
1019  ELSE
1020  a = alx
1021  b = blx
1022  END IF
1023  ALLOCATE (pot0(nr), pot1(nr))
1024  DO ir = 1, nr
1025  IF (rho0(ir) > 1.e-12_dp) THEN
1026  pot0(ir) = 0.5_dp*drho0(ir)/(3._dp*pi*pi*rho0(ir)**fx)
1027  ELSE
1028  pot0(ir) = 0._dp
1029  END IF
1030  END DO
1031  pot1(1:nr) = 1._dp + (a*pot0(1:nr)**2)/(1._dp + b*pot0(1:nr)**2)
1032  pot1(1:nr) = pot1(1:nr)*delta(1:nr)
1033  dfexc = fourpi*integrate_grid(pot1(1:nr), basis0%grid)
1034  DEALLOCATE (pot0, pot1)
1035  CASE DEFAULT
1036  cpabort("Unknown functional in atom_dpot_lda")
1037  END SELECT
1038 
1039  DEALLOCATE (rho0, rho1, drho0, drho1, delta)
1040 
1041  CALL timestop(handle)
1042 
1043  END SUBROUTINE atom_dpot_lda
1044 
1045 ! **************************************************************************************************
1046 !> \brief Initialise object that holds components needed to compute the exchange-correlation
1047 !> functional on the atomic radial grid.
1048 !> \param rho_set electron density object to initialise
1049 !> \param nspins number of spin components
1050 !> \param needs components needed to compute the exchange-correlation functional
1051 !> \param rho electron density on the grid
1052 !> \param drho first derivative of the electron density on the grid
1053 !> \param tau kinetic energy density on the grid
1054 !> \param lap second derivative of the electron density on the grid
1055 !> \param na number of points on the atomic radial grid
1056 !> \par History
1057 !> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
1058 !> * 08.2008 created as calculate_atom_vxc_lda() [Juerg Hutter]
1059 ! **************************************************************************************************
1060  SUBROUTINE fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, na)
1061 
1062  TYPE(xc_rho_set_type), INTENT(INOUT) :: rho_set
1063  INTEGER, INTENT(IN) :: nspins
1064  TYPE(xc_rho_cflags_type), INTENT(in) :: needs
1065  REAL(dp), DIMENSION(:, :), POINTER :: rho, drho, tau, lap
1066  INTEGER, INTENT(IN) :: na
1067 
1068  REAL(kind=dp), PARAMETER :: f13 = (1.0_dp/3.0_dp)
1069 
1070  INTEGER :: ia
1071 
1072  SELECT CASE (nspins)
1073  CASE (1)
1074  cpassert(.NOT. needs%rho_spin)
1075  cpassert(.NOT. needs%drho_spin)
1076  cpassert(.NOT. needs%norm_drho_spin)
1077  cpassert(.NOT. needs%rho_spin_1_3)
1078  cpassert(.NOT. needs%tau_spin)
1079  cpassert(.NOT. needs%drho)
1080  ! Give rho to 1/3
1081  IF (needs%rho_1_3) THEN
1082  DO ia = 1, na
1083  rho_set%rho_1_3(ia, 1, 1) = max(rho(ia, 1), 0.0_dp)**f13
1084  END DO
1085  rho_set%owns%rho_1_3 = .true.
1086  rho_set%has%rho_1_3 = .true.
1087  END IF
1088  ! Give the density
1089  IF (needs%rho) THEN
1090  DO ia = 1, na
1091  rho_set%rho(ia, 1, 1) = rho(ia, 1)
1092  END DO
1093  rho_set%owns%rho = .true.
1094  rho_set%has%rho = .true.
1095  END IF
1096  ! Give the norm of the gradient of the density
1097  IF (needs%norm_drho) THEN
1098  DO ia = 1, na
1099  rho_set%norm_drho(ia, 1, 1) = drho(ia, 1)
1100  END DO
1101  rho_set%owns%norm_drho = .true.
1102  rho_set%has%norm_drho = .true.
1103  END IF
1104  CASE (2)
1105  cpassert(.NOT. needs%drho)
1106  cpassert(.NOT. needs%drho_spin)
1107  ! Give the total density
1108  IF (needs%rho) THEN
1109  DO ia = 1, na
1110  rho_set%rho(ia, 1, 1) = rho(ia, 1) + rho(ia, 2)
1111  END DO
1112  rho_set%owns%rho = .true.
1113  rho_set%has%rho = .true.
1114  END IF
1115  ! Give the norm of the total gradient of the density
1116  IF (needs%norm_drho) THEN
1117  DO ia = 1, na
1118  rho_set%norm_drho(ia, 1, 1) = drho(ia, 1) + drho(ia, 2)
1119  END DO
1120  rho_set%owns%norm_drho = .true.
1121  rho_set%has%norm_drho = .true.
1122  END IF
1123  ! Give rho_spin
1124  IF (needs%rho_spin) THEN
1125  DO ia = 1, na
1126  rho_set%rhoa(ia, 1, 1) = rho(ia, 1)
1127  rho_set%rhob(ia, 1, 1) = rho(ia, 2)
1128  END DO
1129  rho_set%owns%rho_spin = .true.
1130  rho_set%has%rho_spin = .true.
1131  END IF
1132  ! Give rho_spin to 1/3
1133  IF (needs%rho_spin_1_3) THEN
1134  DO ia = 1, na
1135  rho_set%rhoa_1_3(ia, 1, 1) = max(rho(ia, 1), 0.0_dp)**f13
1136  rho_set%rhob_1_3(ia, 1, 1) = max(rho(ia, 2), 0.0_dp)**f13
1137  END DO
1138  rho_set%owns%rho_1_3 = .true.
1139  rho_set%has%rho_1_3 = .true.
1140  END IF
1141  ! Give the norm of the gradient of rhoa and of rhob separatedly
1142  IF (needs%norm_drho_spin) THEN
1143  DO ia = 1, na
1144  rho_set%norm_drhoa(ia, 1, 1) = drho(ia, 1)
1145  rho_set%norm_drhob(ia, 1, 1) = drho(ia, 2)
1146  END DO
1147  rho_set%owns%norm_drho_spin = .true.
1148  rho_set%has%norm_drho_spin = .true.
1149  END IF
1150  !
1151  END SELECT
1152 
1153  ! tau part
1154  IF (needs%tau) THEN
1155  IF (nspins == 2) THEN
1156  DO ia = 1, na
1157  rho_set%tau(ia, 1, 1) = tau(ia, 1) + tau(ia, 2)
1158  END DO
1159  rho_set%owns%tau = .true.
1160  rho_set%has%tau = .true.
1161  ELSE
1162  DO ia = 1, na
1163  rho_set%tau(ia, 1, 1) = tau(ia, 1)
1164  END DO
1165  rho_set%owns%tau = .true.
1166  rho_set%has%tau = .true.
1167  END IF
1168  END IF
1169  IF (needs%tau_spin) THEN
1170  cpassert(nspins == 2)
1171  DO ia = 1, na
1172  rho_set%tau_a(ia, 1, 1) = tau(ia, 1)
1173  rho_set%tau_b(ia, 1, 1) = tau(ia, 2)
1174  END DO
1175  rho_set%owns%tau_spin = .true.
1176  rho_set%has%tau_spin = .true.
1177  END IF
1178 
1179  ! Laplace
1180  IF (needs%laplace_rho) THEN
1181  IF (nspins == 2) THEN
1182  DO ia = 1, na
1183  rho_set%laplace_rho(ia, 1, 1) = lap(ia, 1) + lap(ia, 2)
1184  END DO
1185  rho_set%owns%laplace_rho = .true.
1186  rho_set%has%laplace_rho = .true.
1187  ELSE
1188  DO ia = 1, na
1189  rho_set%laplace_rho(ia, 1, 1) = lap(ia, 1)
1190  END DO
1191  rho_set%owns%laplace_rho = .true.
1192  rho_set%has%laplace_rho = .true.
1193  END IF
1194  END IF
1195  IF (needs%laplace_rho_spin) THEN
1196  cpassert(nspins == 2)
1197  DO ia = 1, na
1198  rho_set%laplace_rhoa(ia, 1, 1) = lap(ia, 1)
1199  rho_set%laplace_rhob(ia, 1, 1) = lap(ia, 2)
1200  END DO
1201  rho_set%owns%laplace_rho_spin = .true.
1202  rho_set%has%laplace_rho_spin = .true.
1203  END IF
1204 
1205  END SUBROUTINE fill_rho_set
1206 
1207 ! **************************************************************************************************
1208 !> \brief Calculate PADE exchange-correlation (xc) energy functional and potential
1209 !> in spin-restricted case.
1210 !> \param rho electron density
1211 !> \param exc XC energy functional
1212 !> \param vxc XC potential
1213 !> \par History
1214 !> * 12.2008 created [Juerg Hutter]
1215 ! **************************************************************************************************
1216  PURE SUBROUTINE lda_pade(rho, exc, vxc)
1217 
1218  REAL(dp), DIMENSION(:), INTENT(IN) :: rho
1219  REAL(dp), DIMENSION(:), INTENT(OUT) :: exc, vxc
1220 
1221  REAL(kind=dp), PARAMETER :: a0 = 0.4581652932831429e+0_dp, a1 = 0.2217058676663745e+1_dp, &
1222  a2 = 0.7405551735357053e+0_dp, a3 = 0.1968227878617998e-1_dp, &
1223  b1 = 1.0000000000000000e+0_dp, b2 = 0.4504130959426697e+1_dp, &
1224  b3 = 0.1110667363742916e+1_dp, b4 = 0.2359291751427506e-1_dp, f13 = (1.0_dp/3.0_dp), &
1225  rsfac = 0.6203504908994000166680065_dp
1226 
1227  INTEGER :: i, n
1228  REAL(kind=dp) :: depade, dpv, dq, epade, p, q, rs
1229 
1230  n = SIZE(rho)
1231  exc(1:n) = 0._dp
1232  vxc(1:n) = 0._dp
1233 
1234  DO i = 1, n
1235  IF (rho(i) > 1.e-20_dp) THEN
1236  rs = rsfac*rho(i)**(-f13)
1237  p = a0 + (a1 + (a2 + a3*rs)*rs)*rs
1238  q = (b1 + (b2 + (b3 + b4*rs)*rs)*rs)*rs
1239  epade = -p/q
1240 
1241  dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs)*rs
1242  dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs)*rs)*rs
1243  depade = f13*rs*(dpv*q - p*dq)/(q*q)
1244 
1245  exc(i) = epade*rho(i)
1246  vxc(i) = epade + depade
1247  END IF
1248  END DO
1249 
1250  END SUBROUTINE lda_pade
1251 
1252 ! **************************************************************************************************
1253 !> \brief Calculate PADE exchange-correlation (xc) energy functional and potential
1254 !> in spin-polarized case.
1255 !> \param rhoa spin-alpha electron density
1256 !> \param rhob spin-beta electron density
1257 !> \param exc XC energy functional
1258 !> \param vxca spin-alpha XC potential
1259 !> \param vxcb spin-beta XC potential
1260 !> \par History
1261 !> * 02.2010 created [Juerg Hutter]
1262 ! **************************************************************************************************
1263  PURE SUBROUTINE lsd_pade(rhoa, rhob, exc, vxca, vxcb)
1264 
1265  REAL(dp), DIMENSION(:), INTENT(IN) :: rhoa, rhob
1266  REAL(dp), DIMENSION(:), INTENT(OUT) :: exc, vxca, vxcb
1267 
1268  REAL(kind=dp), PARAMETER :: a0 = 0.4581652932831429e+0_dp, a1 = 0.2217058676663745e+1_dp, &
1269  a2 = 0.7405551735357053e+0_dp, a3 = 0.1968227878617998e-1_dp, &
1270  b1 = 1.0000000000000000e+0_dp, b2 = 0.4504130959426697e+1_dp, &
1271  b3 = 0.1110667363742916e+1_dp, b4 = 0.2359291751427506e-1_dp, &
1272  da0 = 0.119086804055547e+0_dp, da1 = 0.6157402568883345e+0_dp, &
1273  da2 = 0.1574201515892867e+0_dp, da3 = 0.3532336663397157e-2_dp, &
1274  db1 = 0.0000000000000000e+0_dp, db2 = 0.2673612973836267e+0_dp, &
1275  db3 = 0.2052004607777787e+0_dp, db4 = 0.4200005045691381e-2_dp, f13 = (1.0_dp/3.0_dp), &
1276  f43 = (4.0_dp/3.0_dp)
1277  REAL(kind=dp), PARAMETER :: fxfac = 1.923661050931536319759455_dp, &
1278  rsfac = 0.6203504908994000166680065_dp
1279 
1280  INTEGER :: i, n
1281  REAL(kind=dp) :: dc, dpv, dq, dr, dx, fa0, fa1, fa2, fa3, &
1282  fb1, fb2, fb3, fb4, fx1, fx2, p, q, &
1283  rhoab, rs, x, xp, xq
1284 
1285 ! 1/(2^(4/3) - 2)
1286 
1287  n = SIZE(rhoa)
1288  exc(1:n) = 0._dp
1289  vxca(1:n) = 0._dp
1290  vxcb(1:n) = 0._dp
1291 
1292  DO i = 1, n
1293  rhoab = rhoa(i) + rhob(i)
1294  IF (rhoab > 1.e-20_dp) THEN
1295  rs = rsfac*rhoab**(-f13)
1296 
1297  x = (rhoa(i) - rhob(i))/rhoab
1298  IF (x < -1.0_dp) THEN
1299  fx1 = 1.0_dp
1300  fx2 = -f43*fxfac*2.0_dp**f13
1301  ELSE IF (x > 1.0_dp) THEN
1302  fx1 = 1.0_dp
1303  fx2 = f43*fxfac*2.0_dp**f13
1304  ELSE
1305  fx1 = ((1.0_dp + x)**f43 + (1.0_dp - x)**f43 - 2.0_dp)*fxfac
1306  fx2 = ((1.0_dp + x)**f13 - (1.0_dp - x)**f13)*fxfac*f43
1307  END IF
1308 
1309  fa0 = a0 + fx1*da0
1310  fa1 = a1 + fx1*da1
1311  fa2 = a2 + fx1*da2
1312  fa3 = a3 + fx1*da3
1313  fb1 = b1 + fx1*db1
1314  fb2 = b2 + fx1*db2
1315  fb3 = b3 + fx1*db3
1316  fb4 = b4 + fx1*db4
1317 
1318  p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
1319  q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
1320  dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
1321  dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
1322  4.0_dp*fb4*rs)*rs)*rs
1323  xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
1324  xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
1325 
1326  dr = (dpv*q - p*dq)/(q*q)
1327  dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx2/rhoab
1328  dc = f13*rs*dr - p/q
1329 
1330  exc(i) = -p/q*rhoab
1331  vxca(i) = dc - dx*rhob(i)
1332  vxcb(i) = dc + dx*rhoa(i)
1333  END IF
1334  END DO
1335 
1336  END SUBROUTINE lsd_pade
1337 
1338 END MODULE atom_xc
atom_types
Define the atom type and its sub types.
Definition: atom_types.F:15
atom_utils
Some basic routines for atomic calculations.
Definition: atom_utils.F:15
atom_utils::atom_core_density
subroutine, public atom_core_density(corden, potential, typ, rr)
...
Definition: atom_utils.F:694
atom_utils::atom_density
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
Definition: atom_utils.F:367
atom_utils::numpot_matrix
subroutine, public numpot_matrix(imat, cpot, basis, derivatives)
Calculate a potential matrix by integrating the potential on an atomic radial grid.
Definition: atom_utils.F:1840
atom_utils::coulomb_potential_numeric
subroutine, public coulomb_potential_numeric(cpot, density, grid)
Numerically compute the Coulomb potential on an atomic radial grid.
Definition: atom_utils.F:1077
atom_xc
routines that build the integrals of the Vxc potential calculated for the atomic code
Definition: atom_xc.F:12
atom_xc::atom_dpot_lda
subroutine, public atom_dpot_lda(basis0, pmat0, basis1, pmat1, maxl, functional, dfexc, linxpar)
Estimate the ADMM exchange energy correction term using a model exchange functional.
Definition: atom_xc.F:979
atom_xc::atom_vxc_lsd
subroutine, public atom_vxc_lsd(basis, pmata, pmatb, maxl, xc_section, fexc, xcmata, xcmatb)
Alternative subroutine that calculates atomic exchange-correlation potential in spin-polarized case.
Definition: atom_xc.F:800
atom_xc::calculate_atom_vxc_lda
subroutine, public calculate_atom_vxc_lda(xcmat, atom, xc_section)
Calculate atomic exchange-correlation potential in spin-restricted case.
Definition: atom_xc.F:199
atom_xc::calculate_atom_ext_vxc
subroutine, public calculate_atom_ext_vxc(vxc, atom, lprint, xcmat)
ZMP subroutine for the scf external density from the external v_xc.
Definition: atom_xc.F:148
atom_xc::calculate_atom_zmp
subroutine, public calculate_atom_zmp(ext_density, atom, lprint, xcmat)
ZMP subroutine for the calculation of the effective constraint potential in the atomic code.
Definition: atom_xc.F:78
atom_xc::calculate_atom_vxc_lsd
subroutine, public calculate_atom_vxc_lsd(xcmata, xcmatb, atom, xc_section)
Calculate atomic exchange-correlation potential in spin-polarized case.
Definition: atom_xc.F:404
atom_xc::atom_vxc_lda
subroutine, public atom_vxc_lda(basis, pmat, maxl, xc_section, fexc, xcmat)
Alternative subroutine that calculates atomic exchange-correlation potential in spin-restricted case.
Definition: atom_xc.F:649
atom
Definition: atom.F:9
cp_files
Utility routines to open and close files. Tracking of preconnections.
Definition: cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition: cp_files.F:308
cp_files::get_unit_number
integer function, public get_unit_number(file_name)
Returns the first logical unit that is not preconnected.
Definition: cp_files.F:237
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition: cp_files.F:119
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::xc_none
integer, parameter, public xc_none
Definition: input_constants.F:551
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_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
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
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition: mathconstants.F:110
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition: mathconstants.F:113
xc_atom
Definition: xc_atom.F:9
xc_atom::xc_rho_set_atom_update
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition: xc_atom.F:432
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_tau
integer, parameter, public deriv_tau
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau_b
integer, parameter, public deriv_tau_b
Definition: xc_derivative_desc.F:35
xc_derivative_desc::deriv_tau_a
integer, parameter, public deriv_tau_a
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_zero_all
subroutine, public xc_dset_zero_all(deriv_set)
...
Definition: xc_derivative_set_types.F:167
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_set_types::xc_dset_release
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
Definition: xc_derivative_set_types.F:144
xc_derivative_set_types::xc_dset_create
subroutine, public xc_dset_create(derivative_set, pw_pool, local_bounds)
creates a derivative set object
Definition: xc_derivative_set_types.F:112
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_derivatives
Definition: xc_derivatives.F:9
xc_derivatives::xc_functionals_get_needs
type(xc_rho_cflags_type) function, public xc_functionals_get_needs(functionals, lsd, calc_potential)
...
Definition: xc_derivatives.F:600
xc_derivatives::xc_functionals_eval
subroutine, public xc_functionals_eval(functionals, lsd, rho_set, deriv_set, deriv_order)
...
Definition: xc_derivatives.F:565
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_create
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
Definition: xc_rho_set_types.F:111
xc_rho_set_types::xc_rho_set_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition: xc_rho_set_types.F:129