Line data Source code
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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: deltarho, rho_dum, vxc, vxc1, vxc2
88 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rho
89 :
90 0 : CALL timeset(routineN, handle)
91 :
92 0 : nr = atom%basis%grid%nr
93 0 : z = atom%z
94 0 : ALLOCATE (rho(nr, 1), vxc(nr), vxc1(nr), vxc2(nr), rho_dum(nr), deltarho(nr))
95 0 : CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
96 : !Vxc1
97 0 : int1 = integrate_grid(ext_density, atom%basis%grid)
98 0 : int1 = fourpi*int1
99 0 : CALL coulomb_potential_numeric(vxc1, rho(:, 1), atom%basis%grid)
100 :
101 0 : vxc1 = -vxc1/z
102 :
103 : !Vxc2
104 0 : rho_dum = rho(:, 1)*int1/z
105 0 : deltarho = rho_dum - ext_density
106 0 : int2 = integrate_grid(deltarho, atom%basis%grid)
107 0 : CALL coulomb_potential_numeric(vxc2, deltarho, atom%basis%grid)
108 :
109 0 : vxc2 = vxc2*atom%lambda
110 :
111 : !Vxc
112 0 : vxc = vxc1 + vxc2
113 :
114 0 : atom%energy%exc = fourpi*integrate_grid(vxc, rho(:, 1), atom%basis%grid)
115 0 : atom%rho_diff_integral = fourpi*int2
116 :
117 0 : IF (lprint) THEN
118 0 : extunit = get_unit_number()
119 : CALL open_file(file_name="linear_potentials.dat", file_status="UNKNOWN", &
120 : file_form="FORMATTED", file_action="WRITE", &
121 0 : unit_number=extunit)
122 :
123 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]")')
124 0 : 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 0 : atom%basis%grid%rad(ir), rho(ir, 1), deltarho(ir), vxc(ir), vxc1(ir), vxc2(ir)
127 : END DO
128 0 : CALL close_file(unit_number=extunit)
129 : END IF
130 :
131 0 : IF (PRESENT(xcmat)) CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)
132 :
133 0 : DEALLOCATE (rho, vxc, vxc1, vxc2, rho_dum, deltarho)
134 :
135 0 : CALL timestop(handle)
136 :
137 0 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rho
157 :
158 0 : CALL timeset(routineN, handle)
159 0 : nr = atom%basis%grid%nr
160 :
161 0 : ALLOCATE (rho(nr, 1))
162 :
163 0 : CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
164 :
165 0 : IF (lprint) THEN
166 0 : extunit = get_unit_number()
167 : CALL open_file(file_name="linear_potentials.dat", file_status="UNKNOWN", &
168 : file_form="FORMATTED", file_action="WRITE", &
169 0 : unit_number=extunit)
170 :
171 0 : WRITE (extunit, fmt='("#",T10,"R[bohr]",T36,"Rho[au]",T61,"V_XC[au]")')
172 0 : DO ir = 1, nr
173 : WRITE (extunit, fmt='(T1,E24.15,T26,E24.15,T51,E24.15)') &
174 0 : atom%basis%grid%rad(ir), rho(ir, 1), vxc(ir)
175 : END DO
176 0 : CALL close_file(unit_number=extunit)
177 : END IF
178 :
179 0 : atom%energy%exc = fourpi*integrate_grid(vxc, rho(:, 1), atom%basis%grid)
180 0 : IF (PRESENT(xcmat)) CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)
181 :
182 0 : DEALLOCATE (rho)
183 :
184 0 : CALL timestop(handle)
185 :
186 0 : 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 114490 : 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 114490 : REAL(KIND=dp), DIMENSION(:), POINTER :: exc, vxc
211 114490 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
212 57245 : 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 57245 : CALL timeset(routineN, handle)
220 :
221 57245 : nlcc = .FALSE.
222 57245 : IF (atom%potential%ppot_type == gth_pseudo) THEN
223 50585 : nlcc = atom%potential%gth_pot%nlcc
224 6660 : ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
225 24 : nlcc = atom%potential%upf_pot%core_correction
226 6636 : ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
227 84 : nlcc = atom%potential%sgp_pot%has_nlcc
228 : END IF
229 :
230 57245 : IF (ASSOCIATED(xc_section)) THEN
231 22023 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
232 22023 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
233 :
234 22023 : IF (myfun == xc_none) THEN
235 16 : atom%energy%exc = 0._dp
236 : ELSE
237 22007 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
238 22007 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
239 22007 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
240 :
241 22007 : lsd = .FALSE.
242 22007 : nspins = 1
243 22007 : 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 22007 : nr = atom%basis%grid%nr
251 220070 : bounds(1:2, 1:3) = 1
252 22007 : bounds(2, 1) = nr
253 :
254 : ! create a place where to put the derivatives
255 22007 : 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 22007 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
259 : ! allocate the required 3d arrays where to store rho and drho
260 22007 : CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
261 :
262 22007 : NULLIFY (rho, drho, tau)
263 22007 : IF (needs%rho) THEN
264 66021 : ALLOCATE (rho(nr, 1))
265 22007 : CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
266 22007 : IF (nlcc) THEN
267 437 : CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
268 : END IF
269 : END IF
270 22007 : IF (needs%norm_drho) THEN
271 63783 : ALLOCATE (drho(nr, 1))
272 21261 : CALL atom_density(drho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="DER")
273 21261 : IF (nlcc) THEN
274 437 : CALL atom_core_density(drho(:, 1), atom%potential, typ="DER", rr=atom%basis%grid%rad)
275 : END IF
276 : END IF
277 22007 : IF (needs%tau) THEN
278 324 : ALLOCATE (tau(nr, 1))
279 : CALL atom_density(tau(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
280 108 : typ="KIN", rr=atom%basis%grid%rad2)
281 : END IF
282 22007 : IF (needs%laplace_rho) THEN
283 0 : ALLOCATE (lap(nr, 1))
284 : CALL atom_density(lap(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, &
285 0 : typ="LAP", rr=atom%basis%grid%rad2)
286 0 : IF (nlcc) THEN
287 0 : CALL atom_core_density(lap(:, 1), atom%potential, typ="LAP", rr=atom%basis%grid%rad)
288 : END IF
289 : END IF
290 :
291 22007 : CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
292 :
293 22007 : CALL xc_dset_zero_all(deriv_set)
294 :
295 22007 : deriv_order = 1
296 : CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
297 22007 : deriv_order=deriv_order)
298 :
299 : ! Integration to get the matrix elements and energy
300 22007 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.FALSE.)
301 22007 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
302 22007 : 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 22007 : IF (needs%rho) THEN
314 22007 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], allocate_deriv=.FALSE.)
315 22007 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
316 22007 : 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 22007 : DEALLOCATE (rho)
325 : END IF
326 22007 : IF (needs%norm_drho) THEN
327 21261 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], allocate_deriv=.FALSE.)
328 21261 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
329 21261 : 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 21261 : DEALLOCATE (drho)
338 : END IF
339 22007 : IF (needs%tau) THEN
340 108 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], allocate_deriv=.FALSE.)
341 108 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
342 108 : n1 = SIZE(xcmat%op, 1)
343 108 : n2 = SIZE(xcmat%op, 2)
344 108 : n3 = SIZE(xcmat%op, 3)
345 540 : ALLOCATE (taumat(n1, n2, 0:n3 - 1))
346 919980 : taumat = 0._dp
347 :
348 43308 : xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
349 108 : CALL numpot_matrix(xcmat%op, xcpot(:, 1, 1), atom%basis, 2)
350 43308 : xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
351 108 : CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
352 540 : DO l = 0, 3
353 1226172 : xcmat%op(:, :, l) = xcmat%op(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
354 : END DO
355 :
356 108 : DEALLOCATE (tau)
357 108 : DEALLOCATE (taumat)
358 : END IF
359 :
360 : ! assume lap is not really needed
361 22007 : IF (needs%laplace_rho) THEN
362 0 : DEALLOCATE (lap)
363 : END IF
364 :
365 : ! Release the xc structure used to store the xc derivatives
366 22007 : CALL xc_dset_release(deriv_set)
367 22007 : 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 35222 : nr = atom%basis%grid%nr
375 176110 : ALLOCATE (rho(nr, 1), exc(nr), vxc(nr))
376 :
377 35222 : CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
378 35222 : IF (nlcc) THEN
379 64 : CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
380 : END IF
381 35222 : CALL lda_pade(rho(:, 1), exc, vxc)
382 :
383 35222 : atom%energy%exc = fourpi*integrate_grid(exc, atom%basis%grid)
384 35222 : CALL numpot_matrix(xcmat%op, vxc, atom%basis, 0)
385 :
386 35222 : DEALLOCATE (rho, exc, vxc)
387 :
388 : END IF
389 :
390 57245 : CALL timestop(handle)
391 :
392 1087655 : 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 2670 : 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 2670 : REAL(KIND=dp), DIMENSION(:), POINTER :: exc, vxca, vxcb, xfun
416 2670 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
417 1335 : 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 1335 : CALL timeset(routineN, handle)
425 :
426 1335 : nlcc = .FALSE.
427 1335 : IF (atom%potential%ppot_type == gth_pseudo) THEN
428 1203 : nlcc = atom%potential%gth_pot%nlcc
429 132 : ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
430 0 : nlcc = atom%potential%upf_pot%core_correction
431 0 : CPASSERT(.NOT. nlcc)
432 132 : ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
433 0 : nlcc = atom%potential%sgp_pot%has_nlcc
434 0 : CPASSERT(.NOT. nlcc)
435 : END IF
436 :
437 1335 : IF (ASSOCIATED(xc_section)) THEN
438 921 : IF (nlcc) THEN
439 365 : nr = atom%basis%grid%nr
440 1095 : ALLOCATE (xfun(nr))
441 : END IF
442 :
443 921 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
444 921 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
445 :
446 921 : IF (myfun == xc_none) THEN
447 0 : atom%energy%exc = 0._dp
448 : ELSE
449 921 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
450 921 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
451 921 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
452 :
453 921 : lsd = .TRUE.
454 921 : nspins = 2
455 921 : 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 921 : nr = atom%basis%grid%nr
463 9210 : bounds(1:2, 1:3) = 1
464 921 : bounds(2, 1) = nr
465 :
466 : ! create a place where to put the derivatives
467 921 : 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 921 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
471 : ! allocate the required 3d arrays where to store rho and drho
472 921 : CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
473 :
474 921 : NULLIFY (rho, drho, tau)
475 921 : IF (needs%rho_spin) THEN
476 2763 : ALLOCATE (rho(nr, 2))
477 921 : CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
478 921 : CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
479 921 : IF (nlcc) THEN
480 146365 : xfun = 0.0_dp
481 365 : CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
482 292365 : rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
483 292365 : rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
484 : END IF
485 : END IF
486 921 : IF (needs%norm_drho_spin) THEN
487 2757 : ALLOCATE (drho(nr, 2))
488 919 : CALL atom_density(drho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="DER")
489 919 : CALL atom_density(drho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="DER")
490 919 : IF (nlcc) THEN
491 146365 : xfun = 0.0_dp
492 365 : CALL atom_core_density(xfun(:), atom%potential, typ="DER", rr=atom%basis%grid%rad)
493 292365 : drho(:, 1) = drho(:, 1) + 0.5_dp*xfun(:)
494 292365 : drho(:, 2) = drho(:, 2) + 0.5_dp*xfun(:)
495 : END IF
496 : END IF
497 921 : IF (needs%tau_spin) THEN
498 6 : ALLOCATE (tau(nr, 2))
499 : CALL atom_density(tau(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, &
500 2 : typ="KIN", rr=atom%basis%grid%rad2)
501 : CALL atom_density(tau(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, &
502 2 : typ="KIN", rr=atom%basis%grid%rad2)
503 : END IF
504 921 : IF (needs%laplace_rho_spin) THEN
505 0 : ALLOCATE (lap(nr, 2))
506 : CALL atom_density(lap(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, &
507 0 : typ="LAP", rr=atom%basis%grid%rad2)
508 : CALL atom_density(lap(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, &
509 0 : typ="LAP", rr=atom%basis%grid%rad2)
510 0 : IF (nlcc) THEN
511 0 : xfun = 0.0_dp
512 0 : CALL atom_core_density(xfun(:), atom%potential, typ="LAP", rr=atom%basis%grid%rad)
513 0 : lap(:, 1) = lap(:, 1) + 0.5_dp*xfun(:)
514 0 : lap(:, 2) = lap(:, 2) + 0.5_dp*xfun(:)
515 : END IF
516 : END IF
517 :
518 921 : CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
519 :
520 921 : CALL xc_dset_zero_all(deriv_set)
521 :
522 921 : deriv_order = 1
523 : CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
524 921 : deriv_order=deriv_order)
525 :
526 : ! Integration to get the matrix elements and energy
527 921 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.FALSE.)
528 921 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
529 921 : atom%energy%exc = fourpi*integrate_grid(xcpot(:, 1, 1), atom%basis%grid)
530 :
531 921 : IF (needs%rho_spin) THEN
532 921 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], allocate_deriv=.FALSE.)
533 921 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
534 921 : CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 0)
535 921 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], allocate_deriv=.FALSE.)
536 921 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
537 921 : CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 0)
538 921 : DEALLOCATE (rho)
539 : END IF
540 921 : IF (needs%norm_drho_spin) THEN
541 : ! drhoa
542 919 : NULLIFY (deriv)
543 919 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], allocate_deriv=.FALSE.)
544 919 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
545 919 : CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 1)
546 : ! drhob
547 919 : NULLIFY (deriv)
548 919 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], allocate_deriv=.FALSE.)
549 919 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
550 919 : CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 1)
551 : ! Cross Terms
552 919 : NULLIFY (deriv)
553 919 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
554 919 : IF (ASSOCIATED(deriv)) THEN
555 919 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
556 919 : CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 1)
557 919 : CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 1)
558 : END IF
559 919 : DEALLOCATE (drho)
560 : END IF
561 921 : IF (needs%tau_spin) THEN
562 2 : n1 = SIZE(xcmata%op, 1)
563 2 : n2 = SIZE(xcmata%op, 2)
564 2 : n3 = SIZE(xcmata%op, 3)
565 10 : ALLOCATE (taumat(n1, n2, 0:n3 - 1))
566 :
567 2 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], allocate_deriv=.FALSE.)
568 2 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
569 86 : taumat = 0._dp
570 802 : xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
571 2 : CALL numpot_matrix(xcmata%op, xcpot(:, 1, 1), atom%basis, 2)
572 802 : xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
573 2 : CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
574 10 : DO l = 0, 3
575 106 : xcmata%op(:, :, l) = xcmata%op(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
576 : END DO
577 :
578 2 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], allocate_deriv=.FALSE.)
579 2 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
580 86 : taumat = 0._dp
581 802 : xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
582 2 : CALL numpot_matrix(xcmatb%op, xcpot(:, 1, 1), atom%basis, 2)
583 802 : xcpot(:, 1, 1) = xcpot(:, 1, 1)/atom%basis%grid%rad2(:)
584 2 : CALL numpot_matrix(taumat, xcpot(:, 1, 1), atom%basis, 0)
585 10 : DO l = 0, 3
586 106 : xcmatb%op(:, :, l) = xcmatb%op(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
587 : END DO
588 :
589 2 : DEALLOCATE (tau)
590 2 : DEALLOCATE (taumat)
591 : END IF
592 :
593 : ! assume lap is not really needed
594 921 : IF (needs%laplace_rho_spin) THEN
595 0 : DEALLOCATE (lap)
596 : END IF
597 :
598 : ! Release the xc structure used to store the xc derivatives
599 921 : CALL xc_dset_release(deriv_set)
600 921 : CALL xc_rho_set_release(rho_set)
601 :
602 : END IF !xc_none
603 :
604 921 : IF (nlcc) DEALLOCATE (xfun)
605 :
606 : ELSE
607 :
608 : ! we don't have an xc_section, use a default setup
609 414 : nr = atom%basis%grid%nr
610 2898 : ALLOCATE (rho(nr, 2), exc(nr), vxca(nr), vxcb(nr))
611 414 : IF (nlcc) ALLOCATE (xfun(nr))
612 :
613 414 : CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
614 414 : CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
615 414 : IF (nlcc) THEN
616 0 : xfun(:) = 0.0_dp
617 0 : CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
618 0 : rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
619 0 : rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
620 : END IF
621 414 : CALL lsd_pade(rho(:, 1), rho(:, 2), exc, vxca, vxcb)
622 :
623 414 : atom%energy%exc = fourpi*integrate_grid(exc, atom%basis%grid)
624 414 : CALL numpot_matrix(xcmata%op, vxca, atom%basis, 0)
625 414 : CALL numpot_matrix(xcmatb%op, vxcb, atom%basis, 0)
626 :
627 414 : DEALLOCATE (rho, exc, vxca, vxcb)
628 414 : IF (nlcc) DEALLOCATE (xfun)
629 :
630 : END IF
631 :
632 1335 : CALL timestop(handle)
633 :
634 25365 : 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 8 : 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 16 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
664 8 : 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 8 : CALL timeset(routineN, handle)
672 :
673 8 : CPASSERT(ASSOCIATED(xc_section))
674 :
675 8 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
676 8 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
677 :
678 8 : IF (myfun == xc_none) THEN
679 0 : fexc = 0._dp
680 : ELSE
681 8 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
682 8 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
683 8 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
684 :
685 8 : lsd = .FALSE.
686 8 : nspins = 1
687 8 : 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 8 : nr = basis%grid%nr
695 80 : bounds(1:2, 1:3) = 1
696 8 : bounds(2, 1) = nr
697 :
698 : ! create a place where to put the derivatives
699 8 : 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 8 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
703 : ! allocate the required 3d arrays where to store rho and drho
704 8 : CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
705 :
706 8 : NULLIFY (rho, drho, tau)
707 8 : IF (needs%rho) THEN
708 24 : ALLOCATE (rho(nr, 1))
709 8 : CALL atom_density(rho(:, 1), pmat, basis, maxl, typ="RHO")
710 : END IF
711 8 : IF (needs%norm_drho) THEN
712 24 : ALLOCATE (drho(nr, 1))
713 8 : CALL atom_density(drho(:, 1), pmat, basis, maxl, typ="DER")
714 : END IF
715 8 : IF (needs%tau) THEN
716 0 : ALLOCATE (tau(nr, 1))
717 0 : CALL atom_density(tau(:, 1), pmat, basis, maxl, typ="KIN", rr=basis%grid%rad2)
718 : END IF
719 8 : IF (needs%laplace_rho) THEN
720 0 : ALLOCATE (lap(nr, 1))
721 0 : CALL atom_density(lap(:, 1), pmat, basis, maxl, typ="LAP", rr=basis%grid%rad2)
722 : END IF
723 :
724 8 : CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
725 :
726 8 : CALL xc_dset_zero_all(deriv_set)
727 :
728 8 : deriv_order = 1
729 : CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
730 8 : deriv_order=deriv_order)
731 :
732 : ! Integration to get the matrix elements and energy
733 8 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.FALSE.)
734 8 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
735 8 : fexc = fourpi*integrate_grid(xcpot(:, 1, 1), basis%grid)
736 :
737 8 : IF (needs%rho) THEN
738 8 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rho], allocate_deriv=.FALSE.)
739 8 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
740 8 : CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 0)
741 8 : DEALLOCATE (rho)
742 : END IF
743 8 : IF (needs%norm_drho) THEN
744 8 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho], allocate_deriv=.FALSE.)
745 8 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
746 8 : CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 1)
747 8 : DEALLOCATE (drho)
748 : END IF
749 8 : IF (needs%tau) THEN
750 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau], allocate_deriv=.FALSE.)
751 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
752 0 : n1 = SIZE(xcmat, 1)
753 0 : n2 = SIZE(xcmat, 2)
754 0 : n3 = SIZE(xcmat, 3)
755 0 : ALLOCATE (taumat(n1, n2, 0:n3 - 1))
756 0 : taumat = 0._dp
757 :
758 0 : xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
759 0 : CALL numpot_matrix(xcmat, xcpot(:, 1, 1), basis, 2)
760 0 : xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
761 0 : CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
762 0 : DO l = 0, 3
763 0 : xcmat(:, :, l) = xcmat(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
764 : END DO
765 :
766 0 : DEALLOCATE (tau)
767 0 : DEALLOCATE (taumat)
768 : END IF
769 :
770 : ! we assume lap is not really needed
771 8 : IF (needs%laplace_rho) THEN
772 0 : DEALLOCATE (lap)
773 : END IF
774 :
775 : ! Release the xc structure used to store the xc derivatives
776 8 : CALL xc_dset_release(deriv_set)
777 8 : CALL xc_rho_set_release(rho_set)
778 :
779 : END IF !xc_none
780 :
781 8 : CALL timestop(handle)
782 :
783 152 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: drho, lap, rho, tau
815 0 : 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 0 : CALL timeset(routineN, handle)
823 :
824 0 : CPASSERT(ASSOCIATED(xc_section))
825 :
826 0 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
827 0 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
828 :
829 0 : IF (myfun == xc_none) THEN
830 0 : fexc = 0._dp
831 : ELSE
832 0 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
833 0 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
834 0 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
835 :
836 0 : lsd = .TRUE.
837 0 : nspins = 2
838 0 : 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 0 : nr = basis%grid%nr
846 0 : bounds(1:2, 1:3) = 1
847 0 : bounds(2, 1) = nr
848 :
849 : ! create a place where to put the derivatives
850 0 : 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 0 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
854 : ! allocate the required 3d arrays where to store rho and drho
855 0 : CALL xc_rho_set_atom_update(rho_set, needs, nspins, bounds)
856 :
857 0 : NULLIFY (rho, drho, tau)
858 0 : IF (needs%rho_spin) THEN
859 0 : ALLOCATE (rho(nr, 2))
860 0 : CALL atom_density(rho(:, 1), pmata, basis, maxl, typ="RHO")
861 0 : CALL atom_density(rho(:, 2), pmatb, basis, maxl, typ="RHO")
862 : END IF
863 0 : IF (needs%norm_drho_spin) THEN
864 0 : ALLOCATE (drho(nr, 2))
865 0 : CALL atom_density(drho(:, 1), pmata, basis, maxl, typ="DER")
866 0 : CALL atom_density(drho(:, 2), pmatb, basis, maxl, typ="DER")
867 : END IF
868 0 : IF (needs%tau_spin) THEN
869 0 : ALLOCATE (tau(nr, 2))
870 0 : CALL atom_density(tau(:, 1), pmata, basis, maxl, typ="KIN", rr=basis%grid%rad2)
871 0 : CALL atom_density(tau(:, 2), pmatb, basis, maxl, typ="KIN", rr=basis%grid%rad2)
872 : END IF
873 0 : IF (needs%laplace_rho_spin) THEN
874 0 : ALLOCATE (lap(nr, 2))
875 0 : CALL atom_density(lap(:, 1), pmata, basis, maxl, typ="LAP", rr=basis%grid%rad2)
876 0 : CALL atom_density(lap(:, 2), pmatb, basis, maxl, typ="LAP", rr=basis%grid%rad2)
877 : END IF
878 :
879 0 : CALL fill_rho_set(rho_set, nspins, needs, rho, drho, tau, lap, nr)
880 :
881 0 : CALL xc_dset_zero_all(deriv_set)
882 :
883 0 : deriv_order = 1
884 : CALL xc_functionals_eval(xc_fun_section, lsd=lsd, rho_set=rho_set, deriv_set=deriv_set, &
885 0 : deriv_order=deriv_order)
886 :
887 : ! Integration to get the matrix elements and energy
888 0 : deriv => xc_dset_get_derivative(deriv_set, [INTEGER::], allocate_deriv=.FALSE.)
889 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
890 0 : fexc = fourpi*integrate_grid(xcpot(:, 1, 1), basis%grid)
891 :
892 0 : IF (needs%rho_spin) THEN
893 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhoa], allocate_deriv=.FALSE.)
894 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
895 0 : CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 0)
896 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_rhob], allocate_deriv=.FALSE.)
897 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
898 0 : CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 0)
899 0 : DEALLOCATE (rho)
900 : END IF
901 0 : IF (needs%norm_drho_spin) THEN
902 : ! drhoa
903 0 : NULLIFY (deriv)
904 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhoa], allocate_deriv=.FALSE.)
905 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
906 0 : CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 1)
907 : ! drhob
908 0 : NULLIFY (deriv)
909 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drhob], allocate_deriv=.FALSE.)
910 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
911 0 : CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 1)
912 : ! Cross Terms
913 0 : NULLIFY (deriv)
914 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_norm_drho])
915 0 : IF (ASSOCIATED(deriv)) THEN
916 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
917 0 : CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 1)
918 0 : CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 1)
919 : END IF
920 0 : DEALLOCATE (drho)
921 : END IF
922 0 : IF (needs%tau_spin) THEN
923 0 : n1 = SIZE(xcmata, 1)
924 0 : n2 = SIZE(xcmata, 2)
925 0 : n3 = SIZE(xcmata, 3)
926 0 : ALLOCATE (taumat(n1, n2, 0:n3 - 1))
927 :
928 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_a], allocate_deriv=.FALSE.)
929 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
930 0 : taumat = 0._dp
931 0 : xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
932 0 : CALL numpot_matrix(xcmata, xcpot(:, 1, 1), basis, 2)
933 0 : xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
934 0 : CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
935 0 : DO l = 0, 3
936 0 : xcmata(:, :, l) = xcmata(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
937 : END DO
938 :
939 0 : deriv => xc_dset_get_derivative(deriv_set, [deriv_tau_b], allocate_deriv=.FALSE.)
940 0 : CALL xc_derivative_get(deriv, deriv_data=xcpot)
941 0 : taumat = 0._dp
942 0 : xcpot(:, 1, 1) = 0.5_dp*xcpot(:, 1, 1)
943 0 : CALL numpot_matrix(xcmatb, xcpot(:, 1, 1), basis, 2)
944 0 : xcpot(:, 1, 1) = xcpot(:, 1, 1)/basis%grid%rad2(:)
945 0 : CALL numpot_matrix(taumat, xcpot(:, 1, 1), basis, 0)
946 0 : DO l = 0, 3
947 0 : xcmatb(:, :, l) = xcmatb(:, :, l) + REAL(l*(l + 1), dp)*taumat(:, :, l)
948 : END DO
949 :
950 0 : DEALLOCATE (tau)
951 0 : DEALLOCATE (taumat)
952 : END IF
953 :
954 : ! Release the xc structure used to store the xc derivatives
955 0 : CALL xc_dset_release(deriv_set)
956 0 : CALL xc_rho_set_release(rho_set)
957 :
958 : END IF !xc_none
959 :
960 0 : CALL timestop(handle)
961 :
962 0 : 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 3 : 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 3 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: delta, drho0, drho1, pot0, pot1, rho0, &
995 : rho1
996 :
997 3 : CALL timeset(routineN, handle)
998 :
999 3 : nr = basis0%grid%nr
1000 :
1001 12 : ALLOCATE (rho0(nr), drho0(nr))
1002 3 : CALL atom_density(rho0, pmat0, basis0, maxl, typ="RHO")
1003 3 : CALL atom_density(drho0, pmat0, basis0, maxl, typ="DER")
1004 : !
1005 9 : ALLOCATE (rho1(nr), drho1(nr))
1006 3 : CALL atom_density(rho1, pmat1, basis1, maxl, typ="RHO")
1007 : ! CONSIDER TO REMOVE [?]: drho1 is calculated but it is not used.
1008 3 : CALL atom_density(drho1, pmat1, basis1, maxl, typ="DER")
1009 : !
1010 6 : ALLOCATE (delta(nr))
1011 3 : fx = 4.0_dp/3.0_dp
1012 1203 : delta(1:nr) = fs*(rho0(1:nr)**fx - rho1(1:nr)**fx)
1013 :
1014 6 : SELECT CASE (TRIM(functional))
1015 : CASE ("LINX")
1016 3 : IF (PRESENT(linxpar)) THEN
1017 0 : a = linxpar(1)
1018 0 : b = linxpar(2)
1019 : ELSE
1020 : a = alx
1021 : b = blx
1022 : END IF
1023 9 : ALLOCATE (pot0(nr), pot1(nr))
1024 1203 : DO ir = 1, nr
1025 1203 : IF (rho0(ir) > 1.e-12_dp) THEN
1026 1146 : pot0(ir) = 0.5_dp*drho0(ir)/(3._dp*pi*pi*rho0(ir)**fx)
1027 : ELSE
1028 54 : pot0(ir) = 0._dp
1029 : END IF
1030 : END DO
1031 1203 : pot1(1:nr) = 1._dp + (a*pot0(1:nr)**2)/(1._dp + b*pot0(1:nr)**2)
1032 1203 : pot1(1:nr) = pot1(1:nr)*delta(1:nr)
1033 3 : dfexc = fourpi*integrate_grid(pot1(1:nr), basis0%grid)
1034 3 : DEALLOCATE (pot0, pot1)
1035 : CASE DEFAULT
1036 3 : CPABORT("Unknown functional in atom_dpot_lda")
1037 : END SELECT
1038 :
1039 3 : DEALLOCATE (rho0, rho1, drho0, drho1, delta)
1040 :
1041 3 : CALL timestop(handle)
1042 :
1043 3 : 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 22936 : 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 44951 : SELECT CASE (nspins)
1073 : CASE (1)
1074 22015 : CPASSERT(.NOT. needs%rho_spin)
1075 22015 : CPASSERT(.NOT. needs%drho_spin)
1076 22015 : CPASSERT(.NOT. needs%norm_drho_spin)
1077 22015 : CPASSERT(.NOT. needs%rho_spin_1_3)
1078 22015 : CPASSERT(.NOT. needs%tau_spin)
1079 22015 : CPASSERT(.NOT. needs%drho)
1080 : ! Give rho to 1/3
1081 22015 : IF (needs%rho_1_3) THEN
1082 354095 : DO ia = 1, na
1083 354095 : rho_set%rho_1_3(ia, 1, 1) = MAX(rho(ia, 1), 0.0_dp)**f13
1084 : END DO
1085 895 : rho_set%owns%rho_1_3 = .TRUE.
1086 895 : rho_set%has%rho_1_3 = .TRUE.
1087 : END IF
1088 : ! Give the density
1089 22015 : IF (needs%rho) THEN
1090 9139647 : DO ia = 1, na
1091 9139647 : rho_set%rho(ia, 1, 1) = rho(ia, 1)
1092 : END DO
1093 22015 : rho_set%owns%rho = .TRUE.
1094 22015 : rho_set%has%rho = .TRUE.
1095 : END IF
1096 : ! Give the norm of the gradient of the density
1097 22015 : IF (needs%norm_drho) THEN
1098 8837869 : DO ia = 1, na
1099 8837869 : rho_set%norm_drho(ia, 1, 1) = drho(ia, 1)
1100 : END DO
1101 21269 : rho_set%owns%norm_drho = .TRUE.
1102 21269 : rho_set%has%norm_drho = .TRUE.
1103 : END IF
1104 : CASE (2)
1105 921 : CPASSERT(.NOT. needs%drho)
1106 921 : CPASSERT(.NOT. needs%drho_spin)
1107 : ! Give the total density
1108 921 : IF (needs%rho) THEN
1109 0 : DO ia = 1, na
1110 0 : rho_set%rho(ia, 1, 1) = rho(ia, 1) + rho(ia, 2)
1111 : END DO
1112 0 : rho_set%owns%rho = .TRUE.
1113 0 : rho_set%has%rho = .TRUE.
1114 : END IF
1115 : ! Give the norm of the total gradient of the density
1116 921 : IF (needs%norm_drho) THEN
1117 365719 : DO ia = 1, na
1118 365719 : rho_set%norm_drho(ia, 1, 1) = drho(ia, 1) + drho(ia, 2)
1119 : END DO
1120 919 : rho_set%owns%norm_drho = .TRUE.
1121 919 : rho_set%has%norm_drho = .TRUE.
1122 : END IF
1123 : ! Give rho_spin
1124 921 : IF (needs%rho_spin) THEN
1125 366521 : DO ia = 1, na
1126 365600 : rho_set%rhoa(ia, 1, 1) = rho(ia, 1)
1127 366521 : rho_set%rhob(ia, 1, 1) = rho(ia, 2)
1128 : END DO
1129 921 : rho_set%owns%rho_spin = .TRUE.
1130 921 : rho_set%has%rho_spin = .TRUE.
1131 : END IF
1132 : ! Give rho_spin to 1/3
1133 921 : IF (needs%rho_spin_1_3) THEN
1134 170024 : DO ia = 1, na
1135 169600 : rho_set%rhoa_1_3(ia, 1, 1) = MAX(rho(ia, 1), 0.0_dp)**f13
1136 170024 : rho_set%rhob_1_3(ia, 1, 1) = MAX(rho(ia, 2), 0.0_dp)**f13
1137 : END DO
1138 424 : rho_set%owns%rho_1_3 = .TRUE.
1139 424 : rho_set%has%rho_1_3 = .TRUE.
1140 : END IF
1141 : ! Give the norm of the gradient of rhoa and of rhob separatedly
1142 23857 : IF (needs%norm_drho_spin) THEN
1143 365719 : DO ia = 1, na
1144 364800 : rho_set%norm_drhoa(ia, 1, 1) = drho(ia, 1)
1145 365719 : rho_set%norm_drhob(ia, 1, 1) = drho(ia, 2)
1146 : END DO
1147 919 : rho_set%owns%norm_drho_spin = .TRUE.
1148 919 : rho_set%has%norm_drho_spin = .TRUE.
1149 : END IF
1150 : !
1151 : END SELECT
1152 :
1153 : ! tau part
1154 22936 : IF (needs%tau) THEN
1155 108 : IF (nspins == 2) THEN
1156 0 : DO ia = 1, na
1157 0 : rho_set%tau(ia, 1, 1) = tau(ia, 1) + tau(ia, 2)
1158 : END DO
1159 0 : rho_set%owns%tau = .TRUE.
1160 0 : rho_set%has%tau = .TRUE.
1161 : ELSE
1162 43308 : DO ia = 1, na
1163 43308 : rho_set%tau(ia, 1, 1) = tau(ia, 1)
1164 : END DO
1165 108 : rho_set%owns%tau = .TRUE.
1166 108 : rho_set%has%tau = .TRUE.
1167 : END IF
1168 : END IF
1169 22936 : IF (needs%tau_spin) THEN
1170 2 : CPASSERT(nspins == 2)
1171 802 : DO ia = 1, na
1172 800 : rho_set%tau_a(ia, 1, 1) = tau(ia, 1)
1173 802 : rho_set%tau_b(ia, 1, 1) = tau(ia, 2)
1174 : END DO
1175 2 : rho_set%owns%tau_spin = .TRUE.
1176 2 : rho_set%has%tau_spin = .TRUE.
1177 : END IF
1178 :
1179 : ! Laplace
1180 22936 : IF (needs%laplace_rho) THEN
1181 0 : IF (nspins == 2) THEN
1182 0 : DO ia = 1, na
1183 0 : rho_set%laplace_rho(ia, 1, 1) = lap(ia, 1) + lap(ia, 2)
1184 : END DO
1185 0 : rho_set%owns%laplace_rho = .TRUE.
1186 0 : rho_set%has%laplace_rho = .TRUE.
1187 : ELSE
1188 0 : DO ia = 1, na
1189 0 : rho_set%laplace_rho(ia, 1, 1) = lap(ia, 1)
1190 : END DO
1191 0 : rho_set%owns%laplace_rho = .TRUE.
1192 0 : rho_set%has%laplace_rho = .TRUE.
1193 : END IF
1194 : END IF
1195 22936 : IF (needs%laplace_rho_spin) THEN
1196 0 : CPASSERT(nspins == 2)
1197 0 : DO ia = 1, na
1198 0 : rho_set%laplace_rhoa(ia, 1, 1) = lap(ia, 1)
1199 0 : rho_set%laplace_rhob(ia, 1, 1) = lap(ia, 2)
1200 : END DO
1201 0 : rho_set%owns%laplace_rho_spin = .TRUE.
1202 0 : rho_set%has%laplace_rho_spin = .TRUE.
1203 : END IF
1204 :
1205 22936 : 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 35222 : 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 35222 : n = SIZE(rho)
1231 14130370 : exc(1:n) = 0._dp
1232 14130370 : vxc(1:n) = 0._dp
1233 :
1234 14130370 : DO i = 1, n
1235 14130370 : IF (rho(i) > 1.e-20_dp) THEN
1236 12863668 : rs = rsfac*rho(i)**(-f13)
1237 12863668 : p = a0 + (a1 + (a2 + a3*rs)*rs)*rs
1238 12863668 : q = (b1 + (b2 + (b3 + b4*rs)*rs)*rs)*rs
1239 12863668 : epade = -p/q
1240 :
1241 12863668 : dpv = a1 + (2.0_dp*a2 + 3.0_dp*a3*rs)*rs
1242 12863668 : dq = b1 + (2.0_dp*b2 + (3.0_dp*b3 + 4.0_dp*b4*rs)*rs)*rs
1243 12863668 : depade = f13*rs*(dpv*q - p*dq)/(q*q)
1244 :
1245 12863668 : exc(i) = epade*rho(i)
1246 12863668 : vxc(i) = epade + depade
1247 : END IF
1248 : END DO
1249 :
1250 35222 : 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 414 : 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 414 : n = SIZE(rhoa)
1288 166014 : exc(1:n) = 0._dp
1289 166014 : vxca(1:n) = 0._dp
1290 166014 : vxcb(1:n) = 0._dp
1291 :
1292 166014 : DO i = 1, n
1293 165600 : rhoab = rhoa(i) + rhob(i)
1294 166014 : IF (rhoab > 1.e-20_dp) THEN
1295 149382 : rs = rsfac*rhoab**(-f13)
1296 :
1297 149382 : x = (rhoa(i) - rhob(i))/rhoab
1298 149382 : IF (x < -1.0_dp) THEN
1299 : fx1 = 1.0_dp
1300 : fx2 = -f43*fxfac*2.0_dp**f13
1301 149382 : ELSE IF (x > 1.0_dp) THEN
1302 : fx1 = 1.0_dp
1303 : fx2 = f43*fxfac*2.0_dp**f13
1304 : ELSE
1305 149382 : fx1 = ((1.0_dp + x)**f43 + (1.0_dp - x)**f43 - 2.0_dp)*fxfac
1306 149382 : fx2 = ((1.0_dp + x)**f13 - (1.0_dp - x)**f13)*fxfac*f43
1307 : END IF
1308 :
1309 149382 : fa0 = a0 + fx1*da0
1310 149382 : fa1 = a1 + fx1*da1
1311 149382 : fa2 = a2 + fx1*da2
1312 149382 : fa3 = a3 + fx1*da3
1313 149382 : fb1 = b1 + fx1*db1
1314 149382 : fb2 = b2 + fx1*db2
1315 149382 : fb3 = b3 + fx1*db3
1316 149382 : fb4 = b4 + fx1*db4
1317 :
1318 149382 : p = fa0 + (fa1 + (fa2 + fa3*rs)*rs)*rs
1319 149382 : q = (fb1 + (fb2 + (fb3 + fb4*rs)*rs)*rs)*rs
1320 149382 : dpv = fa1 + (2.0_dp*fa2 + 3.0_dp*fa3*rs)*rs
1321 : dq = fb1 + (2.0_dp*fb2 + (3.0_dp*fb3 + &
1322 149382 : 4.0_dp*fb4*rs)*rs)*rs
1323 149382 : xp = da0 + (da1 + (da2 + da3*rs)*rs)*rs
1324 149382 : xq = (db1 + (db2 + (db3 + db4*rs)*rs)*rs)*rs
1325 :
1326 149382 : dr = (dpv*q - p*dq)/(q*q)
1327 149382 : dx = 2.0_dp*(xp*q - p*xq)/(q*q)*fx2/rhoab
1328 149382 : dc = f13*rs*dr - p/q
1329 :
1330 149382 : exc(i) = -p/q*rhoab
1331 149382 : vxca(i) = dc - dx*rhob(i)
1332 149382 : vxcb(i) = dc + dx*rhoa(i)
1333 : END IF
1334 : END DO
1335 :
1336 414 : END SUBROUTINE lsd_pade
1337 :
1338 : END MODULE atom_xc
|
