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