Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 density in the basis set of spherical primitives
11 : ! **************************************************************************************************
12 : MODULE qs_vxc_atom
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind
15 : USE basis_set_types, ONLY: get_gto_basis_set,&
16 : gto_basis_set_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE input_constants, ONLY: xc_none
19 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
20 : section_vals_type,&
21 : section_vals_val_get
22 : USE kinds, ONLY: dp
23 : USE memory_utilities, ONLY: reallocate
24 : USE message_passing, ONLY: mp_para_env_type
25 : USE orbital_pointers, ONLY: indso,&
26 : nsoset
27 : USE paw_basis_types, ONLY: get_paw_basis_info
28 : USE qs_environment_types, ONLY: get_qs_env,&
29 : qs_environment_type
30 : USE qs_grid_atom, ONLY: grid_atom_type
31 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
32 : harmonics_atom_type
33 : USE qs_kind_types, ONLY: get_qs_kind,&
34 : has_nlcc,&
35 : qs_kind_type
36 : USE qs_linres_types, ONLY: nablavks_atom_type
37 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
38 : rho_atom_coeff,&
39 : rho_atom_type
40 : USE util, ONLY: get_limit
41 : USE xc_atom, ONLY: fill_rho_set,&
42 : vxc_of_r_epr,&
43 : vxc_of_r_new,&
44 : xc_2nd_deriv_of_r,&
45 : xc_rho_set_atom_update
46 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
47 : xc_dset_create,&
48 : xc_dset_release,&
49 : xc_dset_zero_all
50 : USE xc_derivatives, ONLY: xc_functionals_get_needs
51 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
52 : USE xc_rho_set_types, ONLY: xc_rho_set_create,&
53 : xc_rho_set_release,&
54 : xc_rho_set_type
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc_atom'
62 :
63 : PUBLIC :: calculate_vxc_atom, &
64 : calculate_vxc_atom_epr, &
65 : calculate_xc_2nd_deriv_atom, &
66 : calc_rho_angular, &
67 : calculate_gfxc_atom, &
68 : gfxc_atom_diff, &
69 : gaVxcgb_noGC
70 :
71 : CONTAINS
72 :
73 : ! **************************************************************************************************
74 : !> \brief ...
75 : !> \param qs_env ...
76 : !> \param energy_only ...
77 : !> \param exc1 the on-body ex energy contribution
78 : !> \param adiabatic_rescale_factor ...
79 : !> \param kind_set_external provides a non-default kind_set to use
80 : !> \param rho_atom_set_external provides a non-default atomic density set to use
81 : !> \param xc_section_external provides an external non-default XC
82 : ! **************************************************************************************************
83 63768 : SUBROUTINE calculate_vxc_atom(qs_env, energy_only, exc1, &
84 : adiabatic_rescale_factor, kind_set_external, &
85 : rho_atom_set_external, xc_section_external)
86 :
87 : TYPE(qs_environment_type), POINTER :: qs_env
88 : LOGICAL, INTENT(IN) :: energy_only
89 : REAL(dp), INTENT(INOUT) :: exc1
90 : REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
91 : TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
92 : POINTER :: kind_set_external
93 : TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
94 : POINTER :: rho_atom_set_external
95 : TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section_external
96 :
97 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_vxc_atom'
98 :
99 : INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
100 : myfun, na, natom, nr, nspins, num_pe
101 : INTEGER, DIMENSION(2, 3) :: bounds
102 21256 : INTEGER, DIMENSION(:), POINTER :: atom_list
103 : LOGICAL :: donlcc, gradient_f, lsd, nlcc, paw_atom, &
104 : tau_f
105 : REAL(dp) :: density_cut, exc_h, exc_s, gradient_cut, &
106 : my_adiabatic_rescale_factor, tau_cut
107 : REAL(dp), DIMENSION(1, 1, 1) :: tau_d
108 : REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
109 42512 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight
110 42512 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
111 21256 : vtau_s, vxc_h, vxc_s
112 42512 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
113 21256 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
114 : TYPE(dft_control_type), POINTER :: dft_control
115 : TYPE(grid_atom_type), POINTER :: grid_atom
116 : TYPE(gto_basis_set_type), POINTER :: basis_1c
117 : TYPE(harmonics_atom_type), POINTER :: harmonics
118 : TYPE(mp_para_env_type), POINTER :: para_env
119 21256 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
120 21256 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
121 21256 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
122 21256 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom_set
123 : TYPE(rho_atom_type), POINTER :: rho_atom
124 : TYPE(section_vals_type), POINTER :: input, my_xc_section, xc_fun_section
125 : TYPE(xc_derivative_set_type) :: deriv_set
126 : TYPE(xc_rho_cflags_type) :: needs
127 : TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
128 :
129 : ! -------------------------------------------------------------------------
130 :
131 21256 : CALL timeset(routineN, handle)
132 :
133 21256 : NULLIFY (atom_list)
134 21256 : NULLIFY (my_kind_set)
135 21256 : NULLIFY (atomic_kind_set)
136 21256 : NULLIFY (grid_atom)
137 21256 : NULLIFY (harmonics)
138 21256 : NULLIFY (input)
139 21256 : NULLIFY (para_env)
140 21256 : NULLIFY (rho_atom)
141 21256 : NULLIFY (my_rho_atom_set)
142 21256 : NULLIFY (rho_nlcc)
143 :
144 21256 : IF (PRESENT(adiabatic_rescale_factor)) THEN
145 44 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
146 : ELSE
147 21212 : my_adiabatic_rescale_factor = 1.0_dp
148 : END IF
149 :
150 : CALL get_qs_env(qs_env=qs_env, &
151 : dft_control=dft_control, &
152 : para_env=para_env, &
153 : atomic_kind_set=atomic_kind_set, &
154 : qs_kind_set=my_kind_set, &
155 : input=input, &
156 21256 : rho_atom_set=my_rho_atom_set)
157 :
158 21256 : IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
159 21256 : IF (PRESENT(rho_atom_set_external)) my_rho_atom_set => rho_atom_set_external
160 :
161 21256 : nlcc = has_nlcc(my_kind_set)
162 :
163 21256 : my_xc_section => section_vals_get_subs_vals(input, "DFT%XC")
164 :
165 21256 : IF (PRESENT(xc_section_external)) my_xc_section => xc_section_external
166 :
167 21256 : xc_fun_section => section_vals_get_subs_vals(my_xc_section, "XC_FUNCTIONAL")
168 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", &
169 21256 : i_val=myfun)
170 :
171 21256 : IF (myfun == xc_none) THEN
172 3162 : exc1 = 0.0_dp
173 12622 : my_rho_atom_set(:)%exc_h = 0.0_dp
174 12622 : my_rho_atom_set(:)%exc_s = 0.0_dp
175 : ELSE
176 : CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
177 18094 : r_val=density_cut)
178 : CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
179 18094 : r_val=gradient_cut)
180 : CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
181 18094 : r_val=tau_cut)
182 :
183 18094 : lsd = dft_control%lsd
184 18094 : nspins = dft_control%nspins
185 : needs = xc_functionals_get_needs(xc_fun_section, &
186 : lsd=lsd, &
187 18094 : calc_potential=.TRUE.)
188 :
189 18094 : gradient_f = (needs%drho .OR. needs%drho_spin)
190 18094 : tau_f = (needs%tau .OR. needs%tau_spin)
191 :
192 : ! Initialize energy contribution from the one center XC terms to zero
193 18094 : exc1 = 0.0_dp
194 :
195 : ! Nullify some pointers for work-arrays
196 18094 : NULLIFY (rho_h, drho_h, rho_s, drho_s, weight)
197 18094 : NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
198 18094 : NULLIFY (tau_h, tau_s)
199 18094 : NULLIFY (vtau_h, vtau_s)
200 :
201 : ! Here starts the loop over all the atoms
202 :
203 53632 : DO ikind = 1, SIZE(atomic_kind_set)
204 35538 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
205 : CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
206 35538 : harmonics=harmonics, grid_atom=grid_atom)
207 35538 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
208 :
209 35538 : IF (.NOT. paw_atom) CYCLE
210 :
211 32490 : nr = grid_atom%nr
212 32490 : na = grid_atom%ng_sphere
213 :
214 : ! Prepare the structures needed to calculate and store the xc derivatives
215 :
216 : ! Array dimension: here anly one dimensional arrays are used,
217 : ! i.e. only the first column of deriv_data is read.
218 : ! The other to dimensions are set to size equal 1
219 324900 : bounds(1:2, 1:3) = 1
220 32490 : bounds(2, 1) = na
221 32490 : bounds(2, 2) = nr
222 :
223 : ! create a place where to put the derivatives
224 32490 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
225 : ! create the place where to store the argument for the functionals
226 : CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
227 32490 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
228 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
229 32490 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
230 :
231 : ! allocate the required 3d arrays where to store rho and drho
232 32490 : CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
233 32490 : CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
234 :
235 32490 : CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
236 32490 : CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
237 32490 : weight => grid_atom%weight
238 32490 : CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
239 32490 : CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
240 : !
241 32490 : IF (gradient_f) THEN
242 19514 : CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
243 19514 : CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
244 19514 : CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
245 19514 : CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
246 : END IF
247 :
248 32490 : IF (tau_f) THEN
249 492 : CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
250 492 : CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
251 492 : CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
252 492 : CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
253 : END IF
254 :
255 : ! NLCC: prepare rho and drho of the core charge for this KIND
256 32490 : donlcc = .FALSE.
257 32490 : IF (nlcc) THEN
258 52 : NULLIFY (rho_nlcc)
259 52 : rho_nlcc => my_kind_set(ikind)%nlcc_pot
260 52 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
261 : END IF
262 :
263 : ! Distribute the atoms of this kind
264 :
265 32490 : num_pe = para_env%num_pe
266 32490 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
267 :
268 58580 : DO iat = bo(1), bo(2)
269 26090 : iatom = atom_list(iat)
270 :
271 26090 : my_rho_atom_set(iatom)%exc_h = 0.0_dp
272 26090 : my_rho_atom_set(iatom)%exc_s = 0.0_dp
273 :
274 26090 : rho_atom => my_rho_atom_set(iatom)
275 90559783 : rho_h = 0.0_dp
276 90559783 : rho_s = 0.0_dp
277 26090 : IF (gradient_f) THEN
278 15435 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
279 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
280 : rho_rad_s=r_s, drho_rad_h=dr_h, &
281 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
282 15435 : rho_rad_s_d=r_s_d)
283 236499174 : drho_h = 0.0_dp
284 236499174 : drho_s = 0.0_dp
285 : ELSE
286 10655 : NULLIFY (r_h, r_s)
287 10655 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
288 10655 : rho_d = 0.0_dp
289 : END IF
290 26090 : IF (tau_f) THEN
291 : !compute tau on the grid all at once
292 356 : CALL calc_tau_atom(tau_h, tau_s, rho_atom, my_kind_set(ikind), nspins)
293 : ELSE
294 25734 : tau_d = 0.0_dp
295 : END IF
296 :
297 1501130 : DO ir = 1, nr
298 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
299 : ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
300 1475040 : r_h_d, r_s_d, drho_h, drho_s)
301 1501130 : IF (donlcc) THEN
302 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
303 2600 : ir, rho_nlcc(:, 1), rho_h, rho_s, rho_nlcc(:, 2), drho_h, drho_s)
304 : END IF
305 : END DO
306 1501130 : DO ir = 1, nr
307 1501130 : IF (tau_f) THEN
308 19100 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
309 19100 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
310 1455940 : ELSE IF (gradient_f) THEN
311 770090 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
312 770090 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
313 : ELSE
314 685850 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
315 685850 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
316 : END IF
317 : END DO
318 :
319 : !-------------------!
320 : ! hard atom density !
321 : !-------------------!
322 26090 : CALL xc_dset_zero_all(deriv_set)
323 : CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight, &
324 : lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h, energy_only=energy_only, &
325 26090 : adiabatic_rescale_factor=my_adiabatic_rescale_factor)
326 26090 : rho_atom%exc_h = rho_atom%exc_h + exc_h
327 :
328 : !-------------------!
329 : ! soft atom density !
330 : !-------------------!
331 26090 : CALL xc_dset_zero_all(deriv_set)
332 : CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight, &
333 : lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s, energy_only=energy_only, &
334 26090 : adiabatic_rescale_factor=my_adiabatic_rescale_factor)
335 26090 : rho_atom%exc_s = rho_atom%exc_s + exc_s
336 :
337 : ! Add contributions to the exc energy
338 :
339 26090 : exc1 = exc1 + rho_atom%exc_h - rho_atom%exc_s
340 :
341 : ! Integration to get the matrix elements relative to the vxc_atom
342 : ! here the products with the primitives is done: gaVxcgb
343 : ! internal transformation to get the integral in cartesian Gaussians
344 :
345 26090 : IF (.NOT. energy_only) THEN
346 24578 : NULLIFY (int_hh, int_ss)
347 24578 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
348 24578 : IF (gradient_f) THEN
349 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
350 14105 : grid_atom, basis_1c, harmonics, nspins)
351 : ELSE
352 : CALL gaVxcgb_noGC(vxc_h, vxc_s, int_hh, int_ss, &
353 10473 : grid_atom, basis_1c, harmonics, nspins)
354 : END IF
355 24578 : IF (tau_f) THEN
356 : CALL dgaVtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
357 356 : grid_atom, basis_1c, harmonics, nspins)
358 : END IF
359 : END IF ! energy_only
360 58580 : NULLIFY (r_h, r_s, dr_h, dr_s)
361 : END DO ! iat
362 :
363 : ! Release the xc structure used to store the xc derivatives
364 32490 : CALL xc_dset_release(deriv_set)
365 32490 : CALL xc_rho_set_release(rho_set_h)
366 118612 : CALL xc_rho_set_release(rho_set_s)
367 :
368 : END DO ! ikind
369 :
370 18094 : CALL para_env%sum(exc1)
371 :
372 18094 : IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
373 18094 : IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
374 18094 : IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
375 18094 : IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
376 :
377 18094 : IF (gradient_f) THEN
378 11104 : IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
379 11104 : IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
380 11104 : IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
381 11104 : IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
382 : END IF
383 :
384 18094 : IF (tau_f) THEN
385 292 : IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
386 292 : IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
387 292 : IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
388 292 : IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
389 : END IF
390 :
391 : END IF !xc_none
392 :
393 21256 : CALL timestop(handle)
394 :
395 914008 : END SUBROUTINE calculate_vxc_atom
396 :
397 : ! **************************************************************************************************
398 : !> \brief ...
399 : !> \param qs_env ...
400 : !> \param exc1 the on-body ex energy contribution
401 : !> \param gradient_atom_set ...
402 : ! **************************************************************************************************
403 30 : SUBROUTINE calculate_vxc_atom_epr(qs_env, exc1, gradient_atom_set)
404 :
405 : TYPE(qs_environment_type), POINTER :: qs_env
406 : REAL(dp), INTENT(INOUT) :: exc1
407 : TYPE(nablavks_atom_type), DIMENSION(:), POINTER :: gradient_atom_set
408 :
409 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_vxc_atom_epr'
410 :
411 : INTEGER :: bo(2), handle, ia, iat, iatom, idir, &
412 : ikind, ir, ispin, myfun, na, natom, &
413 : nr, nspins, num_pe
414 : INTEGER, DIMENSION(2, 3) :: bounds
415 10 : INTEGER, DIMENSION(:), POINTER :: atom_list
416 : LOGICAL :: donlcc, gradient_f, lsd, nlcc, paw_atom, &
417 : tau_f
418 : REAL(dp) :: density_cut, exc_h, exc_s, gradient_cut, &
419 : tau_cut
420 : REAL(dp), DIMENSION(1, 1, 1) :: tau_d
421 : REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
422 20 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight
423 20 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
424 10 : vtau_s, vxc_h, vxc_s
425 20 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
426 10 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
427 : TYPE(dft_control_type), POINTER :: dft_control
428 : TYPE(grid_atom_type), POINTER :: grid_atom
429 : TYPE(gto_basis_set_type), POINTER :: basis_1c
430 : TYPE(harmonics_atom_type), POINTER :: harmonics
431 : TYPE(mp_para_env_type), POINTER :: para_env
432 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
433 10 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
434 10 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
435 10 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom_set
436 : TYPE(rho_atom_type), POINTER :: rho_atom
437 : TYPE(section_vals_type), POINTER :: input, my_xc_section, xc_fun_section
438 : TYPE(xc_derivative_set_type) :: deriv_set
439 : TYPE(xc_rho_cflags_type) :: needs
440 : TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
441 :
442 : ! -------------------------------------------------------------------------
443 :
444 10 : CALL timeset(routineN, handle)
445 :
446 10 : NULLIFY (atom_list)
447 10 : NULLIFY (my_kind_set)
448 10 : NULLIFY (atomic_kind_set)
449 10 : NULLIFY (grid_atom)
450 10 : NULLIFY (harmonics)
451 10 : NULLIFY (input)
452 10 : NULLIFY (para_env)
453 10 : NULLIFY (rho_atom)
454 10 : NULLIFY (my_rho_atom_set)
455 10 : NULLIFY (rho_nlcc)
456 :
457 : CALL get_qs_env(qs_env=qs_env, &
458 : dft_control=dft_control, &
459 : para_env=para_env, &
460 : atomic_kind_set=atomic_kind_set, &
461 : qs_kind_set=my_kind_set, &
462 : input=input, &
463 10 : rho_atom_set=my_rho_atom_set)
464 :
465 10 : nlcc = has_nlcc(my_kind_set)
466 :
467 : my_xc_section => section_vals_get_subs_vals(input, &
468 10 : "PROPERTIES%LINRES%EPR%PRINT%G_TENSOR%XC")
469 10 : xc_fun_section => section_vals_get_subs_vals(my_xc_section, "XC_FUNCTIONAL")
470 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", &
471 10 : i_val=myfun)
472 :
473 10 : IF (myfun == xc_none) THEN
474 0 : exc1 = 0.0_dp
475 0 : my_rho_atom_set(:)%exc_h = 0.0_dp
476 0 : my_rho_atom_set(:)%exc_s = 0.0_dp
477 : ELSE
478 : CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
479 10 : r_val=density_cut)
480 : CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
481 10 : r_val=gradient_cut)
482 : CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
483 10 : r_val=tau_cut)
484 :
485 10 : lsd = dft_control%lsd
486 10 : nspins = dft_control%nspins
487 : needs = xc_functionals_get_needs(xc_fun_section, &
488 : lsd=lsd, &
489 10 : calc_potential=.TRUE.)
490 :
491 : ! whatever the xc, if epr_xc, drho_spin is needed
492 10 : needs%drho_spin = .TRUE.
493 :
494 10 : gradient_f = (needs%drho .OR. needs%drho_spin)
495 10 : tau_f = (needs%tau .OR. needs%tau_spin)
496 :
497 : ! Initialize energy contribution from the one center XC terms to zero
498 10 : exc1 = 0.0_dp
499 :
500 : ! Nullify some pointers for work-arrays
501 10 : NULLIFY (rho_h, drho_h, rho_s, drho_s, weight)
502 10 : NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
503 10 : NULLIFY (tau_h, tau_s)
504 10 : NULLIFY (vtau_h, vtau_s)
505 :
506 : ! Here starts the loop over all the atoms
507 :
508 30 : DO ikind = 1, SIZE(atomic_kind_set)
509 20 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
510 : CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
511 20 : harmonics=harmonics, grid_atom=grid_atom)
512 20 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
513 :
514 20 : IF (.NOT. paw_atom) CYCLE
515 :
516 20 : nr = grid_atom%nr
517 20 : na = grid_atom%ng_sphere
518 :
519 : ! Prepare the structures needed to calculate and store the xc derivatives
520 :
521 : ! Array dimension: here anly one dimensional arrays are used,
522 : ! i.e. only the first column of deriv_data is read.
523 : ! The other to dimensions are set to size equal 1
524 200 : bounds(1:2, 1:3) = 1
525 20 : bounds(2, 1) = na
526 20 : bounds(2, 2) = nr
527 :
528 : ! create a place where to put the derivatives
529 20 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
530 : ! create the place where to store the argument for the functionals
531 : CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
532 20 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
533 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
534 20 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
535 :
536 : ! allocate the required 3d arrays where to store rho and drho
537 20 : CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
538 20 : CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
539 :
540 20 : CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
541 20 : CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
542 20 : weight => grid_atom%weight
543 20 : CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
544 20 : CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
545 : !
546 : IF (gradient_f) THEN
547 20 : CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
548 20 : CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
549 20 : CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
550 20 : CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
551 : END IF
552 :
553 20 : IF (tau_f) THEN
554 0 : CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
555 0 : CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
556 0 : CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
557 0 : CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
558 : END IF
559 :
560 : ! NLCC: prepare rho and drho of the core charge for this KIND
561 20 : donlcc = .FALSE.
562 20 : IF (nlcc) THEN
563 0 : NULLIFY (rho_nlcc)
564 0 : rho_nlcc => my_kind_set(ikind)%nlcc_pot
565 0 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
566 : END IF
567 :
568 : ! Distribute the atoms of this kind
569 :
570 20 : num_pe = para_env%num_pe
571 20 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
572 :
573 35 : DO iat = bo(1), bo(2)
574 15 : iatom = atom_list(iat)
575 :
576 15 : my_rho_atom_set(iatom)%exc_h = 0.0_dp
577 15 : my_rho_atom_set(iatom)%exc_s = 0.0_dp
578 :
579 15 : rho_atom => my_rho_atom_set(iatom)
580 76545 : rho_h = 0.0_dp
581 76545 : rho_s = 0.0_dp
582 : IF (gradient_f) THEN
583 15 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
584 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
585 : rho_rad_s=r_s, drho_rad_h=dr_h, &
586 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
587 15 : rho_rad_s_d=r_s_d)
588 376545 : drho_h = 0.0_dp
589 376545 : drho_s = 0.0_dp
590 : ELSE
591 : NULLIFY (r_h, r_s)
592 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
593 : rho_d = 0.0_dp
594 : END IF
595 15 : IF (tau_f) THEN
596 : !compute tau on the grid all at once
597 0 : CALL calc_tau_atom(tau_h, tau_s, rho_atom, my_kind_set(ikind), nspins)
598 : ELSE
599 15 : tau_d = 0.0_dp
600 : END IF
601 :
602 765 : DO ir = 1, nr
603 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
604 : ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
605 750 : r_h_d, r_s_d, drho_h, drho_s)
606 765 : IF (donlcc) THEN
607 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
608 0 : ir, rho_nlcc(:, 1), rho_h, rho_s, rho_nlcc(:, 2), drho_h, drho_s)
609 : END IF
610 : END DO
611 765 : DO ir = 1, nr
612 765 : IF (tau_f) THEN
613 0 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
614 0 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
615 : ELSE IF (gradient_f) THEN
616 750 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
617 750 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
618 : ELSE
619 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
620 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
621 : END IF
622 : END DO
623 :
624 : !-------------------!
625 : ! hard atom density !
626 : !-------------------!
627 15 : CALL xc_dset_zero_all(deriv_set)
628 : CALL vxc_of_r_epr(xc_fun_section, rho_set_h, deriv_set, needs, weight, &
629 15 : lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
630 15 : rho_atom%exc_h = rho_atom%exc_h + exc_h
631 :
632 : !-------------------!
633 : ! soft atom density !
634 : !-------------------!
635 15 : CALL xc_dset_zero_all(deriv_set)
636 : CALL vxc_of_r_epr(xc_fun_section, rho_set_s, deriv_set, needs, weight, &
637 15 : lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
638 15 : rho_atom%exc_s = rho_atom%exc_s + exc_s
639 :
640 45 : DO ispin = 1, nspins
641 135 : DO idir = 1, 3
642 4620 : DO ir = 1, nr
643 229590 : DO ia = 1, na
644 : gradient_atom_set(iatom)%nablavks_vec_rad_h(idir, ispin)%r_coef(ir, ia) = &
645 : gradient_atom_set(iatom)%nablavks_vec_rad_h(idir, ispin)%r_coef(ir, ia) &
646 225000 : + vxg_h(idir, ia, ir, ispin)
647 : gradient_atom_set(iatom)%nablavks_vec_rad_s(idir, ispin)%r_coef(ir, ia) = &
648 : gradient_atom_set(iatom)%nablavks_vec_rad_s(idir, ispin)%r_coef(ir, ia) &
649 229500 : + vxg_s(idir, ia, ir, ispin)
650 : END DO ! ia
651 : END DO ! ir
652 : END DO ! idir
653 : END DO ! ispin
654 :
655 : ! Add contributions to the exc energy
656 :
657 15 : exc1 = exc1 + rho_atom%exc_h - rho_atom%exc_s
658 :
659 : ! Integration to get the matrix elements relative to the vxc_atom
660 : ! here the products with the primitives is done: gaVxcgb
661 : ! internal transformation to get the integral in cartesian Gaussians
662 :
663 15 : NULLIFY (int_hh, int_ss)
664 15 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
665 : IF (gradient_f) THEN
666 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
667 15 : grid_atom, basis_1c, harmonics, nspins)
668 : ELSE
669 : CALL gaVxcgb_noGC(vxc_h, vxc_s, int_hh, int_ss, &
670 : grid_atom, basis_1c, harmonics, nspins)
671 : END IF
672 15 : IF (tau_f) THEN
673 : CALL dgaVtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
674 0 : grid_atom, basis_1c, harmonics, nspins)
675 : END IF
676 35 : NULLIFY (r_h, r_s, dr_h, dr_s)
677 : END DO ! iat
678 :
679 : ! Release the xc structure used to store the xc derivatives
680 20 : CALL xc_dset_release(deriv_set)
681 20 : CALL xc_rho_set_release(rho_set_h)
682 70 : CALL xc_rho_set_release(rho_set_s)
683 :
684 : END DO ! ikind
685 :
686 10 : CALL para_env%sum(exc1)
687 :
688 10 : IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
689 10 : IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
690 10 : IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
691 10 : IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
692 :
693 : IF (gradient_f) THEN
694 10 : IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
695 10 : IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
696 10 : IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
697 10 : IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
698 : END IF
699 :
700 10 : IF (tau_f) THEN
701 0 : IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
702 0 : IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
703 0 : IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
704 0 : IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
705 : END IF
706 :
707 : END IF !xc_none
708 :
709 10 : CALL timestop(handle)
710 :
711 430 : END SUBROUTINE calculate_vxc_atom_epr
712 :
713 : ! **************************************************************************************************
714 : !> \brief ...
715 : !> \param rho_atom_set ...
716 : !> \param rho1_atom_set ...
717 : !> \param qs_env ...
718 : !> \param xc_section ...
719 : !> \param para_env ...
720 : !> \param do_tddfpt2 New implementation of TDDFT.
721 : !> \param do_triplet ...
722 : !> \param do_sf ...
723 : !> \param kind_set_external ...
724 : ! **************************************************************************************************
725 20604 : SUBROUTINE calculate_xc_2nd_deriv_atom(rho_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
726 : do_tddfpt2, do_triplet, do_sf, kind_set_external)
727 :
728 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho1_atom_set
729 : TYPE(qs_environment_type), POINTER :: qs_env
730 : TYPE(section_vals_type), POINTER :: xc_section
731 : TYPE(mp_para_env_type), INTENT(IN) :: para_env
732 : LOGICAL, INTENT(IN), OPTIONAL :: do_tddfpt2, do_triplet, do_sf
733 : TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
734 : POINTER :: kind_set_external
735 :
736 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_xc_2nd_deriv_atom'
737 :
738 : INTEGER :: atom, handle, iatom, ikind, ir, na, &
739 : natom, nr, nspins
740 : INTEGER, DIMENSION(2) :: local_loop_limit
741 : INTEGER, DIMENSION(2, 3) :: bounds
742 3434 : INTEGER, DIMENSION(:), POINTER :: atom_list
743 : LOGICAL :: gradient_functional, lsd, my_do_sf, &
744 : paw_atom, scale_rho, tau_f
745 : REAL(KIND=dp) :: density_cut, gradient_cut, rtot, tau_cut
746 : REAL(KIND=dp), CONTIGUOUS, DIMENSION(:, :, :), &
747 3434 : POINTER :: vxc_h, vxc_s
748 : REAL(KIND=dp), DIMENSION(1, 1, 1) :: rtau
749 : REAL(KIND=dp), DIMENSION(1, 1, 1, 1) :: rrho
750 3434 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: weight
751 10302 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho1_h, rho1_s, rho_h, rho_s, tau1_h, &
752 3434 : tau1_s, tau_h, tau_s
753 6868 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: drho1_h, drho1_s, drho_h, drho_s, vxg_h, &
754 3434 : vxg_s
755 3434 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
756 : TYPE(grid_atom_type), POINTER :: grid_atom
757 : TYPE(gto_basis_set_type), POINTER :: basis_1c
758 : TYPE(harmonics_atom_type), POINTER :: harmonics
759 3434 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set, qs_kind_set
760 3434 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr1_h, dr1_s, dr_h, dr_s, int_hh, &
761 3434 : int_ss, r1_h, r1_s, r_h, r_s
762 3434 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r1_h_d, r1_s_d, r_h_d, r_s_d
763 : TYPE(rho_atom_type), POINTER :: rho1_atom, rho_atom
764 : TYPE(section_vals_type), POINTER :: input, xc_fun_section
765 : TYPE(xc_derivative_set_type) :: deriv_set
766 : TYPE(xc_rho_cflags_type) :: needs
767 : TYPE(xc_rho_set_type) :: rho1_set_h, rho1_set_s, rho_set_h, &
768 : rho_set_s
769 :
770 : ! -------------------------------------------------------------------------
771 :
772 3434 : CALL timeset(routineN, handle)
773 :
774 3434 : NULLIFY (qs_kind_set)
775 3434 : NULLIFY (rho_h, rho_s, drho_h, drho_s, weight)
776 3434 : NULLIFY (rho1_h, rho1_s, drho1_h, drho1_s)
777 3434 : NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
778 3434 : NULLIFY (tau_h, tau_s, tau1_h, tau1_s)
779 :
780 : CALL get_qs_env(qs_env=qs_env, &
781 : input=input, &
782 : qs_kind_set=qs_kind_set, &
783 3434 : atomic_kind_set=atomic_kind_set)
784 :
785 3434 : IF (PRESENT(kind_set_external)) THEN
786 398 : my_kind_set => kind_set_external
787 : ELSE
788 3036 : my_kind_set => qs_kind_set
789 : END IF
790 :
791 3434 : CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
792 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", &
793 3434 : r_val=density_cut)
794 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", &
795 3434 : r_val=gradient_cut)
796 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", &
797 3434 : r_val=tau_cut)
798 :
799 3434 : my_do_sf = .FALSE.
800 3434 : IF (PRESENT(do_sf)) my_do_sf = do_sf
801 :
802 : xc_fun_section => section_vals_get_subs_vals(xc_section, &
803 3434 : "XC_FUNCTIONAL")
804 3434 : IF (lsd) THEN
805 124 : nspins = 2
806 : ELSE
807 3310 : nspins = 1
808 : END IF
809 :
810 3434 : scale_rho = .FALSE.
811 3434 : IF (PRESENT(do_tddfpt2) .AND. PRESENT(do_triplet)) THEN
812 1652 : IF (nspins == 1 .AND. do_triplet) THEN
813 182 : lsd = .TRUE.
814 182 : scale_rho = .TRUE.
815 : END IF
816 1782 : ELSEIF (PRESENT(do_triplet)) THEN
817 1578 : IF (nspins == 1 .AND. do_triplet) lsd = .TRUE.
818 : END IF
819 :
820 : needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, &
821 3434 : calc_potential=.TRUE.)
822 3434 : gradient_functional = needs%drho .OR. needs%drho_spin
823 3434 : tau_f = (needs%tau .OR. needs%tau_spin)
824 : IF (tau_f) THEN
825 0 : CPABORT("Tau functionals not implemented for GAPW 2nd derivatives")
826 : ELSE
827 3434 : rtau = 0.0_dp
828 : END IF
829 :
830 : ! Here starts the loop over all the atoms
831 10904 : DO ikind = 1, SIZE(atomic_kind_set)
832 :
833 7470 : NULLIFY (atom_list, harmonics, grid_atom)
834 7470 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
835 : CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
836 7470 : harmonics=harmonics, grid_atom=grid_atom)
837 7470 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
838 7470 : IF (.NOT. paw_atom) CYCLE
839 :
840 7010 : nr = grid_atom%nr
841 7010 : na = grid_atom%ng_sphere
842 :
843 : ! Array dimension: here anly one dimensional arrays are used,
844 : ! i.e. only the first column of deriv_data is read.
845 : ! The other to dimensions are set to size equal 1.
846 70100 : bounds(1:2, 1:3) = 1
847 7010 : bounds(2, 1) = na
848 7010 : bounds(2, 2) = nr
849 :
850 7010 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
851 : CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
852 7010 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
853 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
854 7010 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
855 : CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
856 7010 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
857 : CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
858 7010 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
859 :
860 : ! allocate the required 3d arrays where to store rho and drho
861 7010 : IF (nspins == 1 .AND. .NOT. lsd) THEN
862 6496 : CALL xc_rho_set_atom_update(rho_set_h, needs, 1, bounds)
863 6496 : CALL xc_rho_set_atom_update(rho1_set_h, needs, 1, bounds)
864 6496 : CALL xc_rho_set_atom_update(rho_set_s, needs, 1, bounds)
865 6496 : CALL xc_rho_set_atom_update(rho1_set_s, needs, 1, bounds)
866 : ELSE
867 514 : CALL xc_rho_set_atom_update(rho_set_h, needs, 2, bounds)
868 514 : CALL xc_rho_set_atom_update(rho1_set_h, needs, 2, bounds)
869 514 : CALL xc_rho_set_atom_update(rho_set_s, needs, 2, bounds)
870 514 : CALL xc_rho_set_atom_update(rho1_set_s, needs, 2, bounds)
871 : END IF
872 :
873 : ALLOCATE (rho_h(1:na, 1:nr, 1:nspins), rho1_h(1:na, 1:nr, 1:nspins), &
874 98140 : rho_s(1:na, 1:nr, 1:nspins), rho1_s(1:na, 1:nr, 1:nspins))
875 :
876 49070 : ALLOCATE (vxc_h(1:na, 1:nr, 1:nspins), vxc_s(1:na, 1:nr, 1:nspins))
877 18440536 : vxc_h = 0.0_dp
878 18440536 : vxc_s = 0.0_dp
879 :
880 7010 : weight => grid_atom%weight
881 :
882 7010 : IF (gradient_functional) THEN
883 : ALLOCATE (drho_h(1:4, 1:na, 1:nr, 1:nspins), drho1_h(1:4, 1:na, 1:nr, 1:nspins), &
884 72324 : drho_s(1:4, 1:na, 1:nr, 1:nspins), drho1_s(1:4, 1:na, 1:nr, 1:nspins))
885 41328 : ALLOCATE (vxg_h(1:3, 1:na, 1:nr, 1:nspins), vxg_s(1:3, 1:na, 1:nr, 1:nspins))
886 : ELSE
887 : ALLOCATE (drho_h(1, 1, 1, 1), drho1_h(1, 1, 1, 1), &
888 1844 : drho_s(1, 1, 1, 1), drho1_s(1, 1, 1, 1))
889 1844 : ALLOCATE (vxg_h(1, 1, 1, 1), vxg_s(1, 1, 1, 1))
890 1844 : rrho = 0.0_dp
891 : END IF
892 54108868 : vxg_h = 0.0_dp
893 54108868 : vxg_s = 0.0_dp
894 :
895 : ! parallelization
896 7010 : local_loop_limit = get_limit(natom, para_env%num_pe, para_env%mepos)
897 :
898 11998 : DO iatom = local_loop_limit(1), local_loop_limit(2) !1,natom
899 4988 : atom = atom_list(iatom)
900 :
901 4988 : rho_atom_set(atom)%exc_h = 0.0_dp
902 4988 : rho_atom_set(atom)%exc_s = 0.0_dp
903 4988 : rho1_atom_set(atom)%exc_h = 0.0_dp
904 4988 : rho1_atom_set(atom)%exc_s = 0.0_dp
905 :
906 4988 : rho_atom => rho_atom_set(atom)
907 4988 : rho1_atom => rho1_atom_set(atom)
908 4988 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
909 4988 : NULLIFY (r1_h, r1_s, dr1_h, dr1_s, r1_h_d, r1_s_d)
910 13198760 : rho_h = 0.0_dp
911 13198760 : rho_s = 0.0_dp
912 13198760 : rho1_h = 0.0_dp
913 13198760 : rho1_s = 0.0_dp
914 4988 : IF (gradient_functional) THEN
915 : CALL get_rho_atom(rho_atom=rho_atom, &
916 : rho_rad_h=r_h, rho_rad_s=r_s, &
917 : drho_rad_h=dr_h, drho_rad_s=dr_s, &
918 3692 : rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
919 : CALL get_rho_atom(rho_atom=rho1_atom, &
920 : rho_rad_h=r1_h, rho_rad_s=r1_s, &
921 : drho_rad_h=dr1_h, drho_rad_s=dr1_s, &
922 3692 : rho_rad_h_d=r1_h_d, rho_rad_s_d=r1_s_d)
923 97299044 : drho_h = 0.0_dp; drho_s = 0.0_dp
924 97299044 : drho1_h = 0.0_dp; drho1_s = 0.0_dp
925 : ELSE
926 : CALL get_rho_atom(rho_atom=rho_atom, &
927 1296 : rho_rad_h=r_h, rho_rad_s=r_s)
928 : CALL get_rho_atom(rho_atom=rho1_atom, &
929 1296 : rho_rad_h=r1_h, rho_rad_s=r1_s)
930 : END IF
931 :
932 4988 : rtot = 0.0_dp
933 :
934 254388 : DO ir = 1, nr
935 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_functional, &
936 : ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, &
937 249400 : drho_h, drho_s)
938 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_functional, &
939 : ir, r1_h, r1_s, rho1_h, rho1_s, dr1_h, dr1_s, r1_h_d, r1_s_d, &
940 254388 : drho1_h, drho1_s)
941 : END DO
942 4988 : IF (scale_rho) THEN
943 497640 : rho_h = 2.0_dp*rho_h
944 497640 : rho_s = 2.0_dp*rho_s
945 195 : IF (gradient_functional) THEN
946 1882800 : drho_h = 2.0_dp*drho_h
947 1882800 : drho_s = 2.0_dp*drho_s
948 : END IF
949 : END IF
950 :
951 254388 : DO ir = 1, nr
952 254388 : IF (tau_f) THEN
953 0 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
954 0 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
955 0 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
956 0 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
957 249400 : ELSE IF (gradient_functional) THEN
958 184600 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, rtau, na, ir)
959 184600 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, rtau, na, ir)
960 184600 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, rtau, na, ir)
961 184600 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, rtau, na, ir)
962 : ELSE
963 64800 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rrho, rtau, na, ir)
964 64800 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rrho, rtau, na, ir)
965 64800 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rrho, rtau, na, ir)
966 64800 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rrho, rtau, na, ir)
967 : END IF
968 : END DO
969 :
970 : CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
971 : rho_set=rho_set_h, rho1_set=rho1_set_h, &
972 : deriv_set=deriv_set, &
973 : w=weight, vxc=vxc_h, vxg=vxg_h, do_triplet=do_triplet, &
974 4988 : do_sf=my_do_sf)
975 : CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
976 : rho_set=rho_set_s, rho1_set=rho1_set_s, &
977 : deriv_set=deriv_set, &
978 : w=weight, vxc=vxc_s, vxg=vxg_s, do_triplet=do_triplet, &
979 4988 : do_sf=my_do_sf)
980 :
981 4988 : CALL get_rho_atom(rho_atom=rho1_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
982 4988 : IF (gradient_functional) THEN
983 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
984 3692 : grid_atom, basis_1c, harmonics, nspins)
985 : ELSE
986 : CALL gaVxcgb_noGC(vxc_h, vxc_s, int_hh, int_ss, &
987 1296 : grid_atom, basis_1c, harmonics, nspins)
988 : END IF
989 :
990 11998 : NULLIFY (r_h, r_s, dr_h, dr_s)
991 :
992 : END DO
993 :
994 : ! some cleanup
995 7010 : NULLIFY (weight)
996 7010 : DEALLOCATE (rho_h, rho_s, rho1_h, rho1_s, vxc_h, vxc_s)
997 7010 : DEALLOCATE (drho_h, drho_s, vxg_h, vxg_s)
998 7010 : DEALLOCATE (drho1_h, drho1_s)
999 :
1000 7010 : CALL xc_dset_release(deriv_set)
1001 7010 : CALL xc_rho_set_release(rho_set_h)
1002 7010 : CALL xc_rho_set_release(rho1_set_h)
1003 7010 : CALL xc_rho_set_release(rho_set_s)
1004 25384 : CALL xc_rho_set_release(rho1_set_s)
1005 :
1006 : END DO
1007 :
1008 3434 : CALL timestop(handle)
1009 :
1010 291890 : END SUBROUTINE calculate_xc_2nd_deriv_atom
1011 :
1012 : ! **************************************************************************************************
1013 : !> \brief ...
1014 : !> \param qs_env ...
1015 : !> \param rho0_atom_set ...
1016 : !> \param rho1_atom_set ...
1017 : !> \param rho2_atom_set ...
1018 : !> \param kind_set ...
1019 : !> \param xc_section ...
1020 : !> \param is_triplet ...
1021 : !> \param accuracy ...
1022 : ! **************************************************************************************************
1023 0 : SUBROUTINE calculate_gfxc_atom(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
1024 : kind_set, xc_section, is_triplet, accuracy)
1025 :
1026 : TYPE(qs_environment_type), POINTER :: qs_env
1027 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
1028 : rho2_atom_set
1029 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
1030 : TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
1031 : LOGICAL, INTENT(IN) :: is_triplet
1032 : INTEGER, INTENT(IN) :: accuracy
1033 :
1034 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_gfxc_atom'
1035 : REAL(KIND=dp), PARAMETER :: epsrho = 5.e-4_dp
1036 :
1037 : INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
1038 : istep, mspins, myfun, na, natom, nf, &
1039 : nr, ns, nspins, nstep, num_pe
1040 : INTEGER, DIMENSION(2, 3) :: bounds
1041 0 : INTEGER, DIMENSION(:), POINTER :: atom_list
1042 : LOGICAL :: donlcc, gradient_f, lsd, nlcc, paw_atom, &
1043 : tau_f
1044 : REAL(dp) :: beta, density_cut, exc_h, exc_s, &
1045 : gradient_cut, oeps1, oeps2, tau_cut
1046 : REAL(dp), DIMENSION(1, 1, 1) :: tau_d
1047 : REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
1048 0 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight
1049 0 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
1050 0 : rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
1051 0 : tau_h, tau_s, vtau_h, vtau_s, vxc_h, &
1052 0 : vxc_s
1053 0 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
1054 0 : drho_h, drho_s, vxg_h, vxg_s
1055 : REAL(KIND=dp), DIMENSION(-4:4) :: ak, bl
1056 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1057 : TYPE(dft_control_type), POINTER :: dft_control
1058 : TYPE(grid_atom_type), POINTER :: grid_atom
1059 : TYPE(gto_basis_set_type), POINTER :: basis_1c
1060 : TYPE(harmonics_atom_type), POINTER :: harmonics
1061 : TYPE(mp_para_env_type), POINTER :: para_env
1062 0 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
1063 0 : int_ss, r_h, r_s
1064 0 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
1065 : TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
1066 : TYPE(section_vals_type), POINTER :: xc_fun_section
1067 : TYPE(xc_derivative_set_type) :: deriv_set
1068 : TYPE(xc_rho_cflags_type) :: needs
1069 : TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
1070 :
1071 0 : CALL timeset(routineN, handle)
1072 :
1073 0 : ak = 0.0_dp
1074 0 : bl = 0.0_dp
1075 0 : SELECT CASE (accuracy)
1076 : CASE (:4)
1077 0 : nstep = 2
1078 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
1079 0 : bl(-2:2) = [-1.0_dp, 16.0_dp, -30.0_dp, 16.0_dp, -1.0_dp]/12.0_dp
1080 : CASE (5:7)
1081 0 : nstep = 3
1082 0 : ak(-3:3) = [-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp]/60.0_dp
1083 0 : bl(-3:3) = [2.0_dp, -27.0_dp, 270.0_dp, -490.0_dp, 270.0_dp, -27.0_dp, 2.0_dp]/180.0_dp
1084 : CASE (8:)
1085 0 : nstep = 4
1086 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
1087 0 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
1088 : bl(-4:4) = [-1.0_dp, 128.0_dp/9.0_dp, -112.0_dp, 896.0_dp, -14350.0_dp/9.0_dp, &
1089 0 : 896.0_dp, -112.0_dp, 128.0_dp/9.0_dp, -1.0_dp]/560.0_dp
1090 : END SELECT
1091 0 : oeps1 = 1.0_dp/epsrho
1092 0 : oeps2 = 1.0_dp/(epsrho**2)
1093 :
1094 : CALL get_qs_env(qs_env=qs_env, &
1095 : dft_control=dft_control, &
1096 : para_env=para_env, &
1097 0 : atomic_kind_set=atomic_kind_set)
1098 :
1099 0 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1100 0 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1101 :
1102 0 : IF (myfun == xc_none) THEN
1103 : ! no action needed?
1104 : ELSE
1105 0 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
1106 0 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
1107 0 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
1108 :
1109 0 : nlcc = has_nlcc(kind_set)
1110 0 : lsd = dft_control%lsd
1111 0 : nspins = dft_control%nspins
1112 0 : mspins = nspins
1113 0 : IF (is_triplet) THEN
1114 0 : CPASSERT(nspins == 1)
1115 0 : lsd = .TRUE.
1116 0 : mspins = 2
1117 : END IF
1118 0 : needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.TRUE.)
1119 0 : gradient_f = (needs%drho .OR. needs%drho_spin)
1120 0 : tau_f = (needs%tau .OR. needs%tau_spin)
1121 :
1122 : ! Here starts the loop over all the atoms
1123 0 : DO ikind = 1, SIZE(atomic_kind_set)
1124 0 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1125 : CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
1126 0 : harmonics=harmonics, grid_atom=grid_atom)
1127 0 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1128 :
1129 0 : IF (.NOT. paw_atom) CYCLE
1130 :
1131 0 : nr = grid_atom%nr
1132 0 : na = grid_atom%ng_sphere
1133 :
1134 : ! Prepare the structures needed to calculate and store the xc derivatives
1135 :
1136 : ! Array dimension: here anly one dimensional arrays are used,
1137 : ! i.e. only the first column of deriv_data is read.
1138 : ! The other to dimensions are set to size equal 1
1139 0 : bounds(1:2, 1:3) = 1
1140 0 : bounds(2, 1) = na
1141 0 : bounds(2, 2) = nr
1142 :
1143 : ! create a place where to put the derivatives
1144 0 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
1145 : ! create the place where to store the argument for the functionals
1146 : CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
1147 0 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1148 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
1149 0 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1150 :
1151 : ! allocate the required 3d arrays where to store rho and drho
1152 0 : CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
1153 0 : CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
1154 :
1155 0 : weight => grid_atom%weight
1156 :
1157 : ALLOCATE (rho_h(na, nr, mspins), rho_s(na, nr, mspins), &
1158 : rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
1159 0 : rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
1160 0 : ALLOCATE (vxc_h(na, nr, mspins), vxc_s(na, nr, mspins))
1161 0 : IF (gradient_f) THEN
1162 : ALLOCATE (drho_h(4, na, nr, mspins), drho_s(4, na, nr, mspins), &
1163 : drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
1164 0 : drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
1165 0 : ALLOCATE (vxg_h(3, na, nr, mspins), vxg_s(3, na, nr, mspins))
1166 : END IF
1167 0 : IF (tau_f) THEN
1168 : ALLOCATE (tau_h(na, nr, mspins), tau_s(na, nr, mspins), &
1169 : tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
1170 0 : tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
1171 0 : ALLOCATE (vtau_h(na, nr, mspins), vtau_s(na, nr, mspins))
1172 : END IF
1173 : !
1174 : ! NLCC: prepare rho and drho of the core charge for this KIND
1175 0 : donlcc = .FALSE.
1176 0 : IF (nlcc) THEN
1177 0 : NULLIFY (rho_nlcc)
1178 0 : rho_nlcc => kind_set(ikind)%nlcc_pot
1179 0 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
1180 : END IF
1181 :
1182 : ! Distribute the atoms of this kind
1183 0 : num_pe = para_env%num_pe
1184 0 : bo = get_limit(natom, num_pe, para_env%mepos)
1185 :
1186 0 : DO iat = bo(1), bo(2)
1187 0 : iatom = atom_list(iat)
1188 : !
1189 0 : NULLIFY (int_hh, int_ss)
1190 0 : rho0_atom => rho0_atom_set(iatom)
1191 0 : CALL get_rho_atom(rho_atom=rho0_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1192 0 : ALLOCATE (fint_ss(nspins), fint_hh(nspins))
1193 0 : DO ns = 1, nspins
1194 0 : nf = SIZE(int_ss(ns)%r_coef, 1)
1195 0 : ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
1196 0 : nf = SIZE(int_hh(ns)%r_coef, 1)
1197 0 : ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
1198 : END DO
1199 :
1200 : ! RHO0
1201 0 : rho0_h = 0.0_dp
1202 0 : rho0_s = 0.0_dp
1203 0 : rho0_atom => rho0_atom_set(iatom)
1204 0 : IF (gradient_f) THEN
1205 0 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1206 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1207 0 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1208 0 : drho0_h = 0.0_dp
1209 0 : drho0_s = 0.0_dp
1210 : ELSE
1211 0 : NULLIFY (r_h, r_s)
1212 0 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1213 0 : rho_d = 0.0_dp
1214 : END IF
1215 0 : DO ir = 1, nr
1216 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1217 : ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
1218 0 : r_h_d, r_s_d, drho0_h, drho0_s)
1219 0 : IF (donlcc) THEN
1220 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
1221 0 : ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
1222 : END IF
1223 : END DO
1224 0 : IF (tau_f) THEN
1225 : !compute tau on the grid all at once
1226 0 : CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, kind_set(ikind), nspins)
1227 : ELSE
1228 0 : tau_d = 0.0_dp
1229 : END IF
1230 : ! RHO1
1231 0 : rho1_h = 0.0_dp
1232 0 : rho1_s = 0.0_dp
1233 0 : rho1_atom => rho1_atom_set(iatom)
1234 0 : IF (gradient_f) THEN
1235 0 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1236 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1237 0 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1238 0 : drho1_h = 0.0_dp
1239 0 : drho1_s = 0.0_dp
1240 : ELSE
1241 0 : NULLIFY (r_h, r_s)
1242 0 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1243 : END IF
1244 0 : DO ir = 1, nr
1245 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1246 : ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
1247 0 : r_h_d, r_s_d, drho1_h, drho1_s)
1248 : END DO
1249 0 : IF (tau_f) THEN
1250 : !compute tau on the grid all at once
1251 0 : CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, kind_set(ikind), nspins)
1252 : END IF
1253 : ! RHO2
1254 0 : rho2_atom => rho2_atom_set(iatom)
1255 :
1256 0 : DO istep = -nstep, nstep
1257 :
1258 0 : beta = REAL(istep, KIND=dp)*epsrho
1259 :
1260 0 : IF (is_triplet) THEN
1261 0 : rho_h(:, :, 1) = rho0_h(:, :, 1) + beta*rho1_h(:, :, 1)
1262 0 : rho_h(:, :, 2) = rho0_h(:, :, 1)
1263 0 : rho_h = 0.5_dp*rho_h
1264 0 : rho_s(:, :, 1) = rho0_s(:, :, 1) + beta*rho1_s(:, :, 1)
1265 0 : rho_s(:, :, 2) = rho0_s(:, :, 1)
1266 0 : rho_s = 0.5_dp*rho_s
1267 0 : IF (gradient_f) THEN
1268 0 : drho_h(:, :, :, 1) = drho0_h(:, :, :, 1) + beta*drho1_h(:, :, :, 1)
1269 0 : drho_h(:, :, :, 2) = drho0_h(:, :, :, 1)
1270 0 : drho_h = 0.5_dp*drho_h
1271 0 : drho_s(:, :, :, 1) = drho0_s(:, :, :, 1) + beta*drho1_s(:, :, :, 1)
1272 0 : drho_s(:, :, :, 2) = drho0_s(:, :, :, 1)
1273 0 : drho_s = 0.5_dp*drho_s
1274 : END IF
1275 0 : IF (tau_f) THEN
1276 0 : tau_h(:, :, 1) = tau0_h(:, :, 1) + beta*tau1_h(:, :, 1)
1277 0 : tau_h(:, :, 2) = tau0_h(:, :, 1)
1278 0 : tau_h = 0.5_dp*tau0_h
1279 0 : tau_s(:, :, 1) = tau0_s(:, :, 1) + beta*tau1_s(:, :, 1)
1280 0 : tau_s(:, :, 2) = tau0_s(:, :, 1)
1281 0 : tau_s = 0.5_dp*tau0_s
1282 : END IF
1283 : ELSE
1284 0 : rho_h = rho0_h + beta*rho1_h
1285 0 : rho_s = rho0_s + beta*rho1_s
1286 0 : IF (gradient_f) THEN
1287 0 : drho_h = drho0_h + beta*drho1_h
1288 0 : drho_s = drho0_s + beta*drho1_s
1289 : END IF
1290 0 : IF (tau_f) THEN
1291 0 : tau_h = tau0_h + beta*tau1_h
1292 0 : tau_s = tau0_s + beta*tau1_s
1293 : END IF
1294 : END IF
1295 : !
1296 0 : IF (gradient_f) THEN
1297 : drho_h(4, :, :, :) = SQRT( &
1298 : drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
1299 : drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
1300 0 : drho_h(3, :, :, :)*drho_h(3, :, :, :))
1301 :
1302 : drho_s(4, :, :, :) = SQRT( &
1303 : drho_s(1, :, :, :)*drho_s(1, :, :, :) + &
1304 : drho_s(2, :, :, :)*drho_s(2, :, :, :) + &
1305 0 : drho_s(3, :, :, :)*drho_s(3, :, :, :))
1306 : END IF
1307 :
1308 0 : DO ir = 1, nr
1309 0 : IF (tau_f) THEN
1310 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_h, na, ir)
1311 0 : CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_s, na, ir)
1312 0 : ELSE IF (gradient_f) THEN
1313 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, drho_h, tau_d, na, ir)
1314 0 : CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, drho_s, tau_d, na, ir)
1315 : ELSE
1316 0 : CALL fill_rho_set(rho_set_h, lsd, mspins, needs, rho_h, rho_d, tau_d, na, ir)
1317 0 : CALL fill_rho_set(rho_set_s, lsd, mspins, needs, rho_s, rho_d, tau_d, na, ir)
1318 : END IF
1319 : END DO
1320 :
1321 : ! hard atom density !
1322 0 : CALL xc_dset_zero_all(deriv_set)
1323 : CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight, &
1324 0 : lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
1325 0 : IF (is_triplet) THEN
1326 0 : vxc_h(:, :, 1) = vxc_h(:, :, 1) - vxc_h(:, :, 2)
1327 0 : IF (gradient_f) THEN
1328 0 : vxg_h(:, :, :, 1) = vxg_h(:, :, :, 1) - vxg_h(:, :, :, 2)
1329 : END IF
1330 0 : IF (tau_f) THEN
1331 0 : vtau_h(:, :, 1) = vtau_h(:, :, 1) - vtau_h(:, :, 2)
1332 : END IF
1333 : END IF
1334 : ! soft atom density !
1335 0 : CALL xc_dset_zero_all(deriv_set)
1336 : CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight, &
1337 0 : lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
1338 0 : IF (is_triplet) THEN
1339 0 : vxc_s(:, :, 1) = vxc_s(:, :, 1) - vxc_s(:, :, 2)
1340 0 : IF (gradient_f) THEN
1341 0 : vxg_s(:, :, :, 1) = vxg_s(:, :, :, 1) - vxg_s(:, :, :, 2)
1342 : END IF
1343 0 : IF (tau_f) THEN
1344 0 : vtau_s(:, :, 1) = vtau_s(:, :, 1) - vtau_s(:, :, 2)
1345 : END IF
1346 : END IF
1347 : ! potentials
1348 0 : DO ns = 1, nspins
1349 0 : fint_hh(ns)%r_coef(:, :) = 0.0_dp
1350 0 : fint_ss(ns)%r_coef(:, :) = 0.0_dp
1351 : END DO
1352 0 : IF (gradient_f) THEN
1353 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
1354 0 : grid_atom, basis_1c, harmonics, nspins)
1355 : ELSE
1356 : CALL gaVxcgb_noGC(vxc_h, vxc_s, fint_hh, fint_ss, &
1357 0 : grid_atom, basis_1c, harmonics, nspins)
1358 : END IF
1359 0 : IF (tau_f) THEN
1360 : CALL dgaVtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
1361 0 : grid_atom, basis_1c, harmonics, nspins)
1362 : END IF
1363 : ! first derivative fxc
1364 0 : NULLIFY (int_hh, int_ss)
1365 0 : CALL get_rho_atom(rho_atom=rho1_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1366 0 : DO ns = 1, nspins
1367 0 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
1368 0 : int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
1369 : END DO
1370 : ! second derivative gxc
1371 0 : NULLIFY (int_hh, int_ss)
1372 0 : CALL get_rho_atom(rho_atom=rho2_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1373 0 : DO ns = 1, nspins
1374 0 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_ss(ns)%r_coef(:, :)
1375 0 : int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps2*bl(istep)*fint_hh(ns)%r_coef(:, :)
1376 : END DO
1377 : END DO
1378 : !
1379 0 : DO ns = 1, nspins
1380 0 : DEALLOCATE (fint_ss(ns)%r_coef)
1381 0 : DEALLOCATE (fint_hh(ns)%r_coef)
1382 : END DO
1383 0 : DEALLOCATE (fint_ss, fint_hh)
1384 :
1385 : END DO ! iat
1386 :
1387 : ! Release the xc structure used to store the xc derivatives
1388 0 : CALL xc_dset_release(deriv_set)
1389 0 : CALL xc_rho_set_release(rho_set_h)
1390 0 : CALL xc_rho_set_release(rho_set_s)
1391 :
1392 0 : DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
1393 0 : DEALLOCATE (vxc_h, vxc_s)
1394 0 : IF (gradient_f) THEN
1395 0 : DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
1396 0 : DEALLOCATE (vxg_h, vxg_s)
1397 : END IF
1398 0 : IF (tau_f) THEN
1399 0 : DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
1400 0 : DEALLOCATE (vtau_h, vtau_s)
1401 : END IF
1402 :
1403 : END DO ! ikind
1404 :
1405 : END IF !xc_none
1406 :
1407 0 : CALL timestop(handle)
1408 :
1409 0 : END SUBROUTINE calculate_gfxc_atom
1410 :
1411 : ! **************************************************************************************************
1412 : !> \brief ...
1413 : !> \param qs_env ...
1414 : !> \param rho0_atom_set ...
1415 : !> \param rho1_atom_set ...
1416 : !> \param rho2_atom_set ...
1417 : !> \param kind_set ...
1418 : !> \param xc_section ...
1419 : !> \param is_triplet ...
1420 : !> \param accuracy ...
1421 : !> \param epsrho ...
1422 : ! **************************************************************************************************
1423 156 : SUBROUTINE gfxc_atom_diff(qs_env, rho0_atom_set, rho1_atom_set, rho2_atom_set, &
1424 : kind_set, xc_section, is_triplet, accuracy, epsrho)
1425 :
1426 : TYPE(qs_environment_type), POINTER :: qs_env
1427 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set, &
1428 : rho2_atom_set
1429 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
1430 : TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section
1431 : LOGICAL, INTENT(IN) :: is_triplet
1432 : INTEGER, INTENT(IN) :: accuracy
1433 : REAL(KIND=dp), INTENT(IN) :: epsrho
1434 :
1435 : CHARACTER(LEN=*), PARAMETER :: routineN = 'gfxc_atom_diff'
1436 :
1437 : INTEGER :: bo(2), handle, iat, iatom, ikind, ir, &
1438 : istep, mspins, myfun, na, natom, nf, &
1439 : nr, ns, nspins, nstep, num_pe
1440 : INTEGER, DIMENSION(2, 3) :: bounds
1441 52 : INTEGER, DIMENSION(:), POINTER :: atom_list
1442 : LOGICAL :: donlcc, gradient_f, lsd, nlcc, paw_atom, &
1443 : tau_f
1444 : REAL(dp) :: beta, density_cut, gradient_cut, oeps1, &
1445 : tau_cut
1446 52 : REAL(dp), CONTIGUOUS, DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
1447 : REAL(dp), DIMENSION(1, 1, 1) :: tau_d
1448 : REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
1449 104 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight
1450 104 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho0_h, rho0_s, rho1_h, rho1_s, rho_h, &
1451 104 : rho_s, tau0_h, tau0_s, tau1_h, tau1_s, &
1452 52 : tau_h, tau_s, vtau_h, vtau_s
1453 52 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho0_h, drho0_s, drho1_h, drho1_s, &
1454 104 : drho_h, drho_s, vxg_h, vxg_s
1455 : REAL(KIND=dp), DIMENSION(-4:4) :: ak
1456 52 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1457 : TYPE(dft_control_type), POINTER :: dft_control
1458 : TYPE(grid_atom_type), POINTER :: grid_atom
1459 : TYPE(gto_basis_set_type), POINTER :: basis_1c
1460 : TYPE(harmonics_atom_type), POINTER :: harmonics
1461 : TYPE(mp_para_env_type), POINTER :: para_env
1462 52 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, fint_hh, fint_ss, int_hh, &
1463 52 : int_ss, r_h, r_s
1464 52 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
1465 : TYPE(rho_atom_type), POINTER :: rho0_atom, rho1_atom, rho2_atom
1466 : TYPE(section_vals_type), POINTER :: xc_fun_section
1467 : TYPE(xc_derivative_set_type) :: deriv_set
1468 : TYPE(xc_rho_cflags_type) :: needs
1469 : TYPE(xc_rho_set_type) :: rho1_set_h, rho1_set_s, rho_set_h, &
1470 : rho_set_s
1471 :
1472 52 : CALL timeset(routineN, handle)
1473 :
1474 52 : ak = 0.0_dp
1475 52 : SELECT CASE (accuracy)
1476 : CASE (:4)
1477 0 : nstep = 2
1478 0 : ak(-2:2) = [1.0_dp, -8.0_dp, 0.0_dp, 8.0_dp, -1.0_dp]/12.0_dp
1479 : CASE (5:7)
1480 416 : nstep = 3
1481 416 : ak(-3:3) = [-1.0_dp, 9.0_dp, -45.0_dp, 0.0_dp, 45.0_dp, -9.0_dp, 1.0_dp]/60.0_dp
1482 : CASE (8:)
1483 0 : nstep = 4
1484 : ak(-4:4) = [1.0_dp, -32.0_dp/3.0_dp, 56.0_dp, -224.0_dp, 0.0_dp, &
1485 52 : 224.0_dp, -56.0_dp, 32.0_dp/3.0_dp, -1.0_dp]/280.0_dp
1486 : END SELECT
1487 52 : oeps1 = 1.0_dp/epsrho
1488 :
1489 : CALL get_qs_env(qs_env=qs_env, &
1490 : dft_control=dft_control, &
1491 : para_env=para_env, &
1492 52 : atomic_kind_set=atomic_kind_set)
1493 :
1494 52 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
1495 52 : CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", i_val=myfun)
1496 :
1497 52 : IF (myfun == xc_none) THEN
1498 : ! no action needed?
1499 : ELSE
1500 : ! calculate fxc
1501 : CALL calculate_xc_2nd_deriv_atom(rho0_atom_set, rho1_atom_set, qs_env, xc_section, para_env, &
1502 52 : do_triplet=is_triplet, kind_set_external=kind_set)
1503 :
1504 52 : CALL section_vals_val_get(xc_section, "DENSITY_CUTOFF", r_val=density_cut)
1505 52 : CALL section_vals_val_get(xc_section, "GRADIENT_CUTOFF", r_val=gradient_cut)
1506 52 : CALL section_vals_val_get(xc_section, "TAU_CUTOFF", r_val=tau_cut)
1507 :
1508 52 : nlcc = has_nlcc(kind_set)
1509 52 : lsd = dft_control%lsd
1510 52 : nspins = dft_control%nspins
1511 52 : mspins = nspins
1512 52 : IF (is_triplet) THEN
1513 6 : CPASSERT(nspins == 1)
1514 6 : lsd = .TRUE.
1515 6 : mspins = 2
1516 : END IF
1517 52 : needs = xc_functionals_get_needs(xc_fun_section, lsd=lsd, calc_potential=.TRUE.)
1518 52 : gradient_f = (needs%drho .OR. needs%drho_spin)
1519 52 : tau_f = (needs%tau .OR. needs%tau_spin)
1520 :
1521 : ! Here starts the loop over all the atoms
1522 178 : DO ikind = 1, SIZE(atomic_kind_set)
1523 126 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1524 : CALL get_qs_kind(kind_set(ikind), paw_atom=paw_atom, &
1525 126 : harmonics=harmonics, grid_atom=grid_atom)
1526 126 : CALL get_qs_kind(kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1527 :
1528 126 : IF (.NOT. paw_atom) CYCLE
1529 :
1530 120 : nr = grid_atom%nr
1531 120 : na = grid_atom%ng_sphere
1532 :
1533 : ! Prepare the structures needed to calculate and store the xc derivatives
1534 :
1535 : ! Array dimension: here anly one dimensional arrays are used,
1536 : ! i.e. only the first column of deriv_data is read.
1537 : ! The other to dimensions are set to size equal 1
1538 1200 : bounds(1:2, 1:3) = 1
1539 120 : bounds(2, 1) = na
1540 120 : bounds(2, 2) = nr
1541 :
1542 : ! create a place where to put the derivatives
1543 120 : CALL xc_dset_create(deriv_set, local_bounds=bounds)
1544 : ! create the place where to store the argument for the functionals
1545 : CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
1546 120 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1547 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
1548 120 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1549 : CALL xc_rho_set_create(rho1_set_h, bounds, rho_cutoff=density_cut, &
1550 120 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1551 : CALL xc_rho_set_create(rho1_set_s, bounds, rho_cutoff=density_cut, &
1552 120 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
1553 :
1554 : ! allocate the required 3d arrays where to store rho and drho
1555 120 : CALL xc_rho_set_atom_update(rho_set_h, needs, mspins, bounds)
1556 120 : CALL xc_rho_set_atom_update(rho_set_s, needs, mspins, bounds)
1557 120 : CALL xc_rho_set_atom_update(rho1_set_h, needs, mspins, bounds)
1558 120 : CALL xc_rho_set_atom_update(rho1_set_s, needs, mspins, bounds)
1559 :
1560 120 : weight => grid_atom%weight
1561 :
1562 : ALLOCATE (rho_h(na, nr, nspins), rho_s(na, nr, nspins), &
1563 : rho0_h(na, nr, nspins), rho0_s(na, nr, nspins), &
1564 2400 : rho1_h(na, nr, nspins), rho1_s(na, nr, nspins))
1565 840 : ALLOCATE (vxc_h(na, nr, nspins), vxc_s(na, nr, nspins))
1566 120 : IF (gradient_f) THEN
1567 : ALLOCATE (drho_h(4, na, nr, nspins), drho_s(4, na, nr, nspins), &
1568 : drho0_h(4, na, nr, nspins), drho0_s(4, na, nr, nspins), &
1569 1600 : drho1_h(4, na, nr, nspins), drho1_s(4, na, nr, nspins))
1570 640 : ALLOCATE (vxg_h(3, na, nr, nspins), vxg_s(3, na, nr, nspins))
1571 : END IF
1572 120 : IF (tau_f) THEN
1573 : ALLOCATE (tau_h(na, nr, nspins), tau_s(na, nr, nspins), &
1574 : tau0_h(na, nr, nspins), tau0_s(na, nr, nspins), &
1575 0 : tau1_h(na, nr, nspins), tau1_s(na, nr, nspins))
1576 0 : ALLOCATE (vtau_h(na, nr, nspins), vtau_s(na, nr, nspins))
1577 : END IF
1578 : !
1579 : ! NLCC: prepare rho and drho of the core charge for this KIND
1580 120 : donlcc = .FALSE.
1581 120 : IF (nlcc) THEN
1582 0 : NULLIFY (rho_nlcc)
1583 0 : rho_nlcc => kind_set(ikind)%nlcc_pot
1584 0 : IF (ASSOCIATED(rho_nlcc)) donlcc = .TRUE.
1585 : END IF
1586 :
1587 : ! Distribute the atoms of this kind
1588 120 : num_pe = para_env%num_pe
1589 120 : bo = get_limit(natom, num_pe, para_env%mepos)
1590 :
1591 205 : DO iat = bo(1), bo(2)
1592 85 : iatom = atom_list(iat)
1593 : !
1594 85 : NULLIFY (int_hh, int_ss)
1595 85 : rho0_atom => rho0_atom_set(iatom)
1596 85 : CALL get_rho_atom(rho_atom=rho0_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1597 510 : ALLOCATE (fint_ss(nspins), fint_hh(nspins))
1598 170 : DO ns = 1, nspins
1599 85 : nf = SIZE(int_ss(ns)%r_coef, 1)
1600 340 : ALLOCATE (fint_ss(ns)%r_coef(nf, nf))
1601 85 : nf = SIZE(int_hh(ns)%r_coef, 1)
1602 425 : ALLOCATE (fint_hh(ns)%r_coef(nf, nf))
1603 : END DO
1604 :
1605 : ! RHO0
1606 216920 : rho0_h = 0.0_dp
1607 216920 : rho0_s = 0.0_dp
1608 85 : rho0_atom => rho0_atom_set(iatom)
1609 85 : IF (gradient_f) THEN
1610 57 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1611 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1612 57 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1613 715464 : drho0_h = 0.0_dp
1614 715464 : drho0_s = 0.0_dp
1615 : ELSE
1616 28 : NULLIFY (r_h, r_s)
1617 28 : CALL get_rho_atom(rho_atom=rho0_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1618 28 : rho_d = 0.0_dp
1619 : END IF
1620 4335 : DO ir = 1, nr
1621 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1622 : ir, r_h, r_s, rho0_h, rho0_s, dr_h, dr_s, &
1623 4250 : r_h_d, r_s_d, drho0_h, drho0_s)
1624 4335 : IF (donlcc) THEN
1625 : CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
1626 0 : ir, rho_nlcc(:, 1), rho0_h, rho0_s, rho_nlcc(:, 2), drho0_h, drho0_s)
1627 : END IF
1628 : END DO
1629 85 : IF (tau_f) THEN
1630 : !compute tau on the grid all at once
1631 0 : CALL calc_tau_atom(tau0_h, tau0_s, rho0_atom, kind_set(ikind), nspins)
1632 : ELSE
1633 85 : tau_d = 0.0_dp
1634 : END IF
1635 : ! RHO1
1636 216920 : rho1_h = 0.0_dp
1637 216920 : rho1_s = 0.0_dp
1638 85 : rho1_atom => rho1_atom_set(iatom)
1639 85 : IF (gradient_f) THEN
1640 57 : NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1641 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
1642 57 : drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1643 715464 : drho1_h = 0.0_dp
1644 715464 : drho1_s = 0.0_dp
1645 : ELSE
1646 28 : NULLIFY (r_h, r_s)
1647 28 : CALL get_rho_atom(rho_atom=rho1_atom, rho_rad_h=r_h, rho_rad_s=r_s)
1648 : END IF
1649 4335 : DO ir = 1, nr
1650 : CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
1651 : ir, r_h, r_s, rho1_h, rho1_s, dr_h, dr_s, &
1652 4335 : r_h_d, r_s_d, drho1_h, drho1_s)
1653 : END DO
1654 85 : IF (tau_f) THEN
1655 : !compute tau on the grid all at once
1656 0 : CALL calc_tau_atom(tau1_h, tau1_s, rho1_atom, kind_set(ikind), nspins)
1657 : END IF
1658 :
1659 4335 : DO ir = 1, nr
1660 4335 : IF (tau_f) THEN
1661 0 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau1_h, na, ir)
1662 0 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau1_s, na, ir)
1663 4250 : ELSE IF (gradient_f) THEN
1664 2850 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, drho1_h, tau_d, na, ir)
1665 2850 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, drho1_s, tau_d, na, ir)
1666 : ELSE
1667 1400 : CALL fill_rho_set(rho1_set_h, lsd, nspins, needs, rho1_h, rho_d, tau_d, na, ir)
1668 1400 : CALL fill_rho_set(rho1_set_s, lsd, nspins, needs, rho1_s, rho_d, tau_d, na, ir)
1669 : END IF
1670 : END DO
1671 :
1672 : ! RHO2
1673 85 : rho2_atom => rho2_atom_set(iatom)
1674 :
1675 680 : DO istep = -nstep, nstep
1676 :
1677 595 : beta = REAL(istep, KIND=dp)*epsrho
1678 :
1679 3036285 : rho_h = rho0_h + beta*rho1_h
1680 3036285 : rho_s = rho0_s + beta*rho1_s
1681 595 : IF (gradient_f) THEN
1682 10016097 : drho_h = drho0_h + beta*drho1_h
1683 10016097 : drho_s = drho0_s + beta*drho1_s
1684 : END IF
1685 595 : IF (tau_f) THEN
1686 0 : tau_h = tau0_h + beta*tau1_h
1687 0 : tau_s = tau0_s + beta*tau1_s
1688 : END IF
1689 : !
1690 595 : IF (gradient_f) THEN
1691 : drho_h(4, :, :, :) = SQRT( &
1692 : drho_h(1, :, :, :)*drho_h(1, :, :, :) + &
1693 : drho_h(2, :, :, :)*drho_h(2, :, :, :) + &
1694 1018248 : drho_h(3, :, :, :)*drho_h(3, :, :, :))
1695 :
1696 : drho_s(4, :, :, :) = SQRT( &
1697 : drho_s(1, :, :, :)*drho_s(1, :, :, :) + &
1698 : drho_s(2, :, :, :)*drho_s(2, :, :, :) + &
1699 1018248 : drho_s(3, :, :, :)*drho_s(3, :, :, :))
1700 : END IF
1701 :
1702 30345 : DO ir = 1, nr
1703 30345 : IF (tau_f) THEN
1704 0 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
1705 0 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
1706 29750 : ELSE IF (gradient_f) THEN
1707 19950 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
1708 19950 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
1709 : ELSE
1710 9800 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
1711 9800 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
1712 : END IF
1713 : END DO
1714 :
1715 : ! hard atom density !
1716 595 : CALL xc_dset_zero_all(deriv_set)
1717 : CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
1718 : rho_set=rho_set_h, rho1_set=rho1_set_h, &
1719 : deriv_set=deriv_set, &
1720 595 : w=weight, vxc=vxc_h, vxg=vxg_h, do_triplet=is_triplet)
1721 : ! soft atom density !
1722 595 : CALL xc_dset_zero_all(deriv_set)
1723 : CALL xc_2nd_deriv_of_r(xc_section=xc_section, &
1724 : rho_set=rho_set_s, rho1_set=rho1_set_s, &
1725 : deriv_set=deriv_set, &
1726 595 : w=weight, vxc=vxc_s, vxg=vxg_s, do_triplet=is_triplet)
1727 : ! potentials
1728 1190 : DO ns = 1, nspins
1729 1184855 : fint_hh(ns)%r_coef(:, :) = 0.0_dp
1730 1185450 : fint_ss(ns)%r_coef(:, :) = 0.0_dp
1731 : END DO
1732 595 : IF (gradient_f) THEN
1733 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, fint_hh, fint_ss, &
1734 399 : grid_atom, basis_1c, harmonics, nspins)
1735 : ELSE
1736 : CALL gaVxcgb_noGC(vxc_h, vxc_s, fint_hh, fint_ss, &
1737 196 : grid_atom, basis_1c, harmonics, nspins)
1738 : END IF
1739 595 : IF (tau_f) THEN
1740 : CALL dgaVtaudgb(vtau_h, vtau_s, fint_hh, fint_ss, &
1741 0 : grid_atom, basis_1c, harmonics, nspins)
1742 : END IF
1743 : ! second derivative gxc
1744 595 : NULLIFY (int_hh, int_ss)
1745 595 : CALL get_rho_atom(rho_atom=rho2_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
1746 1275 : DO ns = 1, nspins
1747 2369115 : int_ss(ns)%r_coef(:, :) = int_ss(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_ss(ns)%r_coef(:, :)
1748 2369710 : int_hh(ns)%r_coef(:, :) = int_hh(ns)%r_coef(:, :) + oeps1*ak(istep)*fint_hh(ns)%r_coef(:, :)
1749 : END DO
1750 : END DO
1751 : !
1752 170 : DO ns = 1, nspins
1753 85 : DEALLOCATE (fint_ss(ns)%r_coef)
1754 170 : DEALLOCATE (fint_hh(ns)%r_coef)
1755 : END DO
1756 205 : DEALLOCATE (fint_ss, fint_hh)
1757 :
1758 : END DO ! iat
1759 :
1760 : ! Release the xc structure used to store the xc derivatives
1761 120 : CALL xc_dset_release(deriv_set)
1762 120 : CALL xc_rho_set_release(rho_set_h)
1763 120 : CALL xc_rho_set_release(rho_set_s)
1764 120 : CALL xc_rho_set_release(rho1_set_h)
1765 120 : CALL xc_rho_set_release(rho1_set_s)
1766 :
1767 120 : DEALLOCATE (rho_h, rho_s, rho0_h, rho0_s, rho1_h, rho1_s)
1768 120 : DEALLOCATE (vxc_h, vxc_s)
1769 120 : IF (gradient_f) THEN
1770 80 : DEALLOCATE (drho_h, drho_s, drho0_h, drho0_s, drho1_h, drho1_s)
1771 80 : DEALLOCATE (vxg_h, vxg_s)
1772 : END IF
1773 418 : IF (tau_f) THEN
1774 0 : DEALLOCATE (tau_h, tau_s, tau0_h, tau0_s, tau1_h, tau1_s)
1775 0 : DEALLOCATE (vtau_h, vtau_s)
1776 : END IF
1777 :
1778 : END DO ! ikind
1779 :
1780 : END IF !xc_none
1781 :
1782 52 : CALL timestop(handle)
1783 :
1784 3796 : END SUBROUTINE gfxc_atom_diff
1785 :
1786 : ! **************************************************************************************************
1787 : !> \brief ...
1788 : !> \param grid_atom ...
1789 : !> \param harmonics ...
1790 : !> \param nspins ...
1791 : !> \param grad_func ...
1792 : !> \param ir ...
1793 : !> \param r_h ...
1794 : !> \param r_s ...
1795 : !> \param rho_h ...
1796 : !> \param rho_s ...
1797 : !> \param dr_h ...
1798 : !> \param dr_s ...
1799 : !> \param r_h_d ...
1800 : !> \param r_s_d ...
1801 : !> \param drho_h ...
1802 : !> \param drho_s ...
1803 : ! **************************************************************************************************
1804 1983090 : SUBROUTINE calc_rho_angular(grid_atom, harmonics, nspins, grad_func, &
1805 : ir, r_h, r_s, rho_h, rho_s, &
1806 : dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
1807 :
1808 : TYPE(grid_atom_type), POINTER :: grid_atom
1809 : TYPE(harmonics_atom_type), POINTER :: harmonics
1810 : INTEGER, INTENT(IN) :: nspins
1811 : LOGICAL, INTENT(IN) :: grad_func
1812 : INTEGER, INTENT(IN) :: ir
1813 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: r_h, r_s
1814 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
1815 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s
1816 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
1817 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
1818 :
1819 : INTEGER :: ia, iso, ispin, na
1820 : REAL(KIND=dp) :: rad, urad
1821 :
1822 1983090 : CPASSERT(ASSOCIATED(r_h))
1823 1983090 : CPASSERT(ASSOCIATED(r_s))
1824 1983090 : CPASSERT(ASSOCIATED(rho_h))
1825 1983090 : CPASSERT(ASSOCIATED(rho_s))
1826 1983090 : IF (grad_func) THEN
1827 1164840 : CPASSERT(ASSOCIATED(dr_h))
1828 1164840 : CPASSERT(ASSOCIATED(dr_s))
1829 1164840 : CPASSERT(ASSOCIATED(r_h_d))
1830 1164840 : CPASSERT(ASSOCIATED(r_s_d))
1831 1164840 : CPASSERT(ASSOCIATED(drho_h))
1832 1164840 : CPASSERT(ASSOCIATED(drho_s))
1833 : END IF
1834 :
1835 1983090 : na = grid_atom%ng_sphere
1836 1983090 : rad = grid_atom%rad(ir)
1837 1983090 : urad = grid_atom%oorad2l(ir, 1)
1838 4283320 : DO ispin = 1, nspins
1839 34074990 : DO iso = 1, harmonics%max_iso_not0
1840 1523216920 : DO ia = 1, na
1841 : rho_h(ia, ir, ispin) = rho_h(ia, ir, ispin) + &
1842 1491125020 : r_h(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
1843 : rho_s(ia, ir, ispin) = rho_s(ia, ir, ispin) + &
1844 1520916690 : r_s(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
1845 : END DO ! ia
1846 : END DO ! iso
1847 : END DO ! ispin
1848 :
1849 1983090 : IF (grad_func) THEN
1850 2500170 : DO ispin = 1, nspins
1851 20246540 : DO iso = 1, harmonics%max_iso_not0
1852 907941720 : DO ia = 1, na
1853 :
1854 : ! components of the gradient of rho1 hard
1855 : drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + &
1856 : dr_h(ispin)%r_coef(ir, iso)* &
1857 : harmonics%a(1, ia)*harmonics%slm(ia, iso) + &
1858 : r_h_d(1, ispin)%r_coef(ir, iso)* &
1859 888860020 : harmonics%slm(ia, iso)
1860 :
1861 : drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + &
1862 : dr_h(ispin)%r_coef(ir, iso)* &
1863 : harmonics%a(2, ia)*harmonics%slm(ia, iso) + &
1864 : r_h_d(2, ispin)%r_coef(ir, iso)* &
1865 888860020 : harmonics%slm(ia, iso)
1866 :
1867 : drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + &
1868 : dr_h(ispin)%r_coef(ir, iso)* &
1869 : harmonics%a(3, ia)*harmonics%slm(ia, iso) + &
1870 : r_h_d(3, ispin)%r_coef(ir, iso)* &
1871 888860020 : harmonics%slm(ia, iso)
1872 :
1873 : ! components of the gradient of rho1 soft
1874 : drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + &
1875 : dr_s(ispin)%r_coef(ir, iso)* &
1876 : harmonics%a(1, ia)*harmonics%slm(ia, iso) + &
1877 : r_s_d(1, ispin)%r_coef(ir, iso)* &
1878 888860020 : harmonics%slm(ia, iso)
1879 :
1880 : drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + &
1881 : dr_s(ispin)%r_coef(ir, iso)* &
1882 : harmonics%a(2, ia)*harmonics%slm(ia, iso) + &
1883 : r_s_d(2, ispin)%r_coef(ir, iso)* &
1884 888860020 : harmonics%slm(ia, iso)
1885 :
1886 : drho_s(3, ia, ir, ispin) = drho_s(3, ia, ir, ispin) + &
1887 : dr_s(ispin)%r_coef(ir, iso)* &
1888 : harmonics%a(3, ia)*harmonics%slm(ia, iso) + &
1889 : r_s_d(3, ispin)%r_coef(ir, iso)* &
1890 888860020 : harmonics%slm(ia, iso)
1891 :
1892 : drho_h(4, ia, ir, ispin) = SQRT( &
1893 : drho_h(1, ia, ir, ispin)*drho_h(1, ia, ir, ispin) + &
1894 : drho_h(2, ia, ir, ispin)*drho_h(2, ia, ir, ispin) + &
1895 888860020 : drho_h(3, ia, ir, ispin)*drho_h(3, ia, ir, ispin))
1896 :
1897 : drho_s(4, ia, ir, ispin) = SQRT( &
1898 : drho_s(1, ia, ir, ispin)*drho_s(1, ia, ir, ispin) + &
1899 : drho_s(2, ia, ir, ispin)*drho_s(2, ia, ir, ispin) + &
1900 906606390 : drho_s(3, ia, ir, ispin)*drho_s(3, ia, ir, ispin))
1901 :
1902 : END DO ! ia
1903 : END DO ! iso
1904 : END DO ! ispin
1905 : END IF
1906 :
1907 1983090 : END SUBROUTINE calc_rho_angular
1908 :
1909 : ! **************************************************************************************************
1910 : !> \brief Computes tau hard and soft on the atomic grids for meta-GGA calculations
1911 : !> \param tau_h the hard part of tau
1912 : !> \param tau_s the soft part of tau
1913 : !> \param rho_atom ...
1914 : !> \param qs_kind ...
1915 : !> \param nspins ...
1916 : !> \note This is a rewrite to correct a meta-GGA GAPW bug. This is more brute force than the original,
1917 : !> which was done along in qs_rho_atom_methods.F, but makes sure that no corner is cut in
1918 : !> terms of accuracy (A. Bussy)
1919 : ! **************************************************************************************************
1920 356 : SUBROUTINE calc_tau_atom(tau_h, tau_s, rho_atom, qs_kind, nspins)
1921 :
1922 : REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: tau_h, tau_s
1923 : TYPE(rho_atom_type), POINTER :: rho_atom
1924 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
1925 : INTEGER, INTENT(IN) :: nspins
1926 :
1927 : CHARACTER(len=*), PARAMETER :: routineN = 'calc_tau_atom'
1928 :
1929 : INTEGER :: dir, handle, ia, ip, ipgf, ir, iset, &
1930 : iso, ispin, jp, jpgf, jset, jso, l, &
1931 : maxso, na, nr, nset, starti, startj
1932 356 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, o2nindex
1933 : REAL(dp) :: cpc_h, cpc_s
1934 356 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fga, fgr, r1, r2
1935 356 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: a1, a2
1936 356 : REAL(dp), DIMENSION(:, :), POINTER :: slm, zet
1937 356 : REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
1938 : TYPE(grid_atom_type), POINTER :: grid_atom
1939 : TYPE(gto_basis_set_type), POINTER :: basis_1c
1940 : TYPE(harmonics_atom_type), POINTER :: harmonics
1941 :
1942 356 : NULLIFY (harmonics, grid_atom, basis_1c, lmax, lmin, npgf, zet, slm, dslm_dxyz, o2nindex)
1943 :
1944 356 : CALL timeset(routineN, handle)
1945 :
1946 : !We need to put 0.5* grad_g1 dot grad_gw on the grid
1947 : !For this we need grid info, basis info, and projector info
1948 356 : CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
1949 356 : CALL get_qs_kind(qs_kind, basis_set=basis_1c, basis_type="GAPW_1C")
1950 :
1951 356 : nr = grid_atom%nr
1952 356 : na = grid_atom%ng_sphere
1953 :
1954 356 : slm => harmonics%slm
1955 356 : dslm_dxyz => harmonics%dslm_dxyz
1956 :
1957 356 : CALL get_paw_basis_info(basis_1c, o2nindex=o2nindex)
1958 :
1959 : !zeroing tau, assuming it is already allocated
1960 1068412 : tau_h = 0.0_dp
1961 1068412 : tau_s = 0.0_dp
1962 :
1963 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, npgf=npgf, &
1964 356 : nset=nset, zet=zet, maxso=maxso)
1965 :
1966 : !Separate the functions into purely r and purely angular parts, precompute them all
1967 2492 : ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3))
1968 2136 : ALLOCATE (r1(nr, nset*maxso), r2(nr, nset*maxso))
1969 4521212 : a1 = 0.0_dp; a2 = 0.0_dp
1970 1521156 : r1 = 0.0_dp; r2 = 0.0_dp
1971 :
1972 1168 : DO iset = 1, nset
1973 3632 : DO ipgf = 1, npgf(iset)
1974 2464 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1975 10307 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
1976 7031 : l = indso(1, iso)
1977 :
1978 : ! The x derivative of the spherical orbital, divided in angular and radial parts
1979 : ! Two of each are needed because d/dx(r^l Y_lm) * exp(-al*r^2) + r^l Y_lm * ! d/dx(exp-al*r^2)
1980 :
1981 : ! the purely radial part of d/dx(r^l Y_lm) * exp(-al*r^2) (same for y, z)
1982 376781 : r1(1:nr, starti + iso) = grid_atom%rad(1:nr)**(l - 1)*EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
1983 :
1984 : ! the purely radial part of r^l Y_lm * d/dx(exp-al*r^2) (same for y, z)
1985 : r2(1:nr, starti + iso) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1) &
1986 376781 : *EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
1987 :
1988 30588 : DO dir = 1, 3
1989 : ! the purely angular part of d/dx(r^l Y_lm) * exp(-al*r^2)
1990 1115055 : a1(1:na, starti + iso, dir) = dslm_dxyz(dir, 1:na, iso)
1991 :
1992 : ! the purely angular part of r^l Y_lm * d/dx(exp-al*r^2)
1993 1122086 : a2(1:na, starti + iso, dir) = harmonics%a(dir, 1:na)*slm(1:na, iso)
1994 : END DO
1995 :
1996 : END DO !iso
1997 : END DO !ipgf
1998 : END DO !iset
1999 :
2000 : !Compute the matrix products
2001 1068 : ALLOCATE (fga(na, 1))
2002 1068 : ALLOCATE (fgr(nr, 1))
2003 38904 : fga = 0.0_dp; fgr = 0.0_dp
2004 :
2005 1168 : DO iset = 1, nset
2006 3576 : DO jset = 1, nset
2007 10518 : DO ipgf = 1, npgf(iset)
2008 7298 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
2009 36446 : DO jpgf = 1, npgf(jset)
2010 26740 : startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
2011 :
2012 110858 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
2013 376637 : DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
2014 :
2015 273077 : ip = o2nindex(starti + iso)
2016 273077 : jp = o2nindex(startj + jso)
2017 :
2018 : !Two component per function => 4 terms in total
2019 :
2020 : ! take r1*a1(dir) * r1*a1(dir)
2021 14245427 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r1(1:nr, startj + jso)
2022 1092308 : DO dir = 1, 3
2023 42468741 : fga(1:na, 1) = a1(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
2024 :
2025 1911539 : DO ispin = 1, nspins
2026 : !get the projectors
2027 819231 : cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
2028 819231 : cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
2029 :
2030 : !compute contribution to tau
2031 43555512 : DO ir = 1, nr
2032 2207384781 : DO ia = 1, na
2033 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
2034 2164648500 : fgr(ir, 1)*fga(ia, 1)
2035 :
2036 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
2037 2206565550 : fgr(ir, 1)*fga(ia, 1)
2038 : END DO
2039 : END DO
2040 :
2041 : END DO !ispin
2042 : END DO !dir
2043 :
2044 : ! add r1*a1(dir) * r2*a2(dir)
2045 14245427 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r2(1:nr, startj + jso)
2046 1092308 : DO dir = 1, 3
2047 42468741 : fga(1:na, 1) = a1(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
2048 :
2049 1911539 : DO ispin = 1, nspins
2050 : !get the projectors
2051 819231 : cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
2052 819231 : cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
2053 :
2054 : !compute contribution to tau
2055 43555512 : DO ir = 1, nr
2056 2207384781 : DO ia = 1, na
2057 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
2058 2164648500 : fgr(ir, 1)*fga(ia, 1)
2059 :
2060 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
2061 2206565550 : fgr(ir, 1)*fga(ia, 1)
2062 : END DO
2063 : END DO
2064 :
2065 : END DO !ispin
2066 : END DO !dir
2067 :
2068 : ! add r2*a2(dir) * V * r1*a1(dir)
2069 14245427 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r1(1:nr, startj + jso)
2070 1092308 : DO dir = 1, 3
2071 42468741 : fga(1:na, 1) = a2(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
2072 :
2073 1911539 : DO ispin = 1, nspins
2074 : !get the projectors
2075 819231 : cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
2076 819231 : cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
2077 :
2078 : !compute contribution to tau
2079 43555512 : DO ir = 1, nr
2080 2207384781 : DO ia = 1, na
2081 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
2082 2164648500 : fgr(ir, 1)*fga(ia, 1)
2083 :
2084 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
2085 2206565550 : fgr(ir, 1)*fga(ia, 1)
2086 : END DO
2087 : END DO
2088 :
2089 : END DO !ispin
2090 : END DO !dir
2091 :
2092 : ! add the last term: r2*a2(dir) * r2*a2(dir)
2093 14245427 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r2(1:nr, startj + jso)
2094 1169128 : DO dir = 1, 3
2095 42468741 : fga(1:na, 1) = a2(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
2096 :
2097 1911539 : DO ispin = 1, nspins
2098 : !get the projectors
2099 819231 : cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
2100 819231 : cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
2101 :
2102 : !compute contribution to tau
2103 43555512 : DO ir = 1, nr
2104 2207384781 : DO ia = 1, na
2105 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
2106 2164648500 : fgr(ir, 1)*fga(ia, 1)
2107 :
2108 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
2109 2206565550 : fgr(ir, 1)*fga(ia, 1)
2110 : END DO
2111 : END DO
2112 :
2113 : END DO !ispin
2114 : END DO !dir
2115 :
2116 : END DO !jso
2117 : END DO !iso
2118 :
2119 : END DO !jpgf
2120 : END DO !ipgf
2121 : END DO !jset
2122 : END DO !iset
2123 :
2124 356 : DEALLOCATE (o2nindex)
2125 :
2126 356 : CALL timestop(handle)
2127 :
2128 1068 : END SUBROUTINE calc_tau_atom
2129 :
2130 : ! **************************************************************************************************
2131 : !> \brief ...
2132 : !> \param grid_atom ...
2133 : !> \param nspins ...
2134 : !> \param grad_func ...
2135 : !> \param ir ...
2136 : !> \param rho_nlcc ...
2137 : !> \param rho_h ...
2138 : !> \param rho_s ...
2139 : !> \param drho_nlcc ...
2140 : !> \param drho_h ...
2141 : !> \param drho_s ...
2142 : ! **************************************************************************************************
2143 2600 : SUBROUTINE calc_rho_nlcc(grid_atom, nspins, grad_func, &
2144 2600 : ir, rho_nlcc, rho_h, rho_s, drho_nlcc, drho_h, drho_s)
2145 :
2146 : TYPE(grid_atom_type), POINTER :: grid_atom
2147 : INTEGER, INTENT(IN) :: nspins
2148 : LOGICAL, INTENT(IN) :: grad_func
2149 : INTEGER, INTENT(IN) :: ir
2150 : REAL(KIND=dp), DIMENSION(:) :: rho_nlcc
2151 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
2152 : REAL(KIND=dp), DIMENSION(:) :: drho_nlcc
2153 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
2154 :
2155 : INTEGER :: ia, ispin, na
2156 : REAL(KIND=dp) :: drho, dx, dy, dz, rad, rho, urad, xsp
2157 :
2158 2600 : CPASSERT(ASSOCIATED(rho_h))
2159 2600 : CPASSERT(ASSOCIATED(rho_s))
2160 2600 : IF (grad_func) THEN
2161 2600 : CPASSERT(ASSOCIATED(drho_h))
2162 2600 : CPASSERT(ASSOCIATED(drho_s))
2163 : END IF
2164 :
2165 2600 : na = grid_atom%ng_sphere
2166 2600 : rad = grid_atom%rad(ir)
2167 2600 : urad = grid_atom%oorad2l(ir, 1)
2168 :
2169 2600 : xsp = REAL(nspins, KIND=dp)
2170 2600 : rho = rho_nlcc(ir)/xsp
2171 5200 : DO ispin = 1, nspins
2172 132600 : rho_h(1:na, ir, ispin) = rho_h(1:na, ir, ispin) + rho
2173 135200 : rho_s(1:na, ir, ispin) = rho_s(1:na, ir, ispin) + rho
2174 : END DO ! ispin
2175 :
2176 2600 : IF (grad_func) THEN
2177 2600 : drho = drho_nlcc(ir)/xsp
2178 5200 : DO ispin = 1, nspins
2179 135200 : DO ia = 1, na
2180 130000 : IF (grid_atom%azi(ia) == 0.0_dp) THEN
2181 : dx = 0.0_dp
2182 : dy = 0.0_dp
2183 : ELSE
2184 117000 : dx = grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
2185 117000 : dy = grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
2186 : END IF
2187 130000 : dz = grid_atom%cos_pol(ia)
2188 : ! components of the gradient of rho1 hard
2189 130000 : drho_h(1, ia, ir, ispin) = drho_h(1, ia, ir, ispin) + drho*dx
2190 130000 : drho_h(2, ia, ir, ispin) = drho_h(2, ia, ir, ispin) + drho*dy
2191 130000 : drho_h(3, ia, ir, ispin) = drho_h(3, ia, ir, ispin) + drho*dz
2192 : ! components of the gradient of rho1 soft
2193 130000 : drho_s(1, ia, ir, ispin) = drho_s(1, ia, ir, ispin) + drho*dx
2194 130000 : drho_s(2, ia, ir, ispin) = drho_s(2, ia, ir, ispin) + drho*dy
2195 130000 : drho_s(3, ia, ir, ispin) = drho_s(3, ia, ir, ispin) + drho*dz
2196 : ! norm of gradient
2197 : drho_h(4, ia, ir, ispin) = SQRT( &
2198 : drho_h(1, ia, ir, ispin)*drho_h(1, ia, ir, ispin) + &
2199 : drho_h(2, ia, ir, ispin)*drho_h(2, ia, ir, ispin) + &
2200 130000 : drho_h(3, ia, ir, ispin)*drho_h(3, ia, ir, ispin))
2201 :
2202 : drho_s(4, ia, ir, ispin) = SQRT( &
2203 : drho_s(1, ia, ir, ispin)*drho_s(1, ia, ir, ispin) + &
2204 : drho_s(2, ia, ir, ispin)*drho_s(2, ia, ir, ispin) + &
2205 132600 : drho_s(3, ia, ir, ispin)*drho_s(3, ia, ir, ispin))
2206 : END DO ! ia
2207 : END DO ! ispin
2208 : END IF
2209 :
2210 2600 : END SUBROUTINE calc_rho_nlcc
2211 :
2212 : ! **************************************************************************************************
2213 : !> \brief ...
2214 : !> \param vxc_h ...
2215 : !> \param vxc_s ...
2216 : !> \param int_hh ...
2217 : !> \param int_ss ...
2218 : !> \param grid_atom ...
2219 : !> \param basis_1c ...
2220 : !> \param harmonics ...
2221 : !> \param nspins ...
2222 : ! **************************************************************************************************
2223 11965 : SUBROUTINE gaVxcgb_noGC(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
2224 :
2225 : REAL(dp), DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
2226 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
2227 : TYPE(grid_atom_type), POINTER :: grid_atom
2228 : TYPE(gto_basis_set_type), POINTER :: basis_1c
2229 : TYPE(harmonics_atom_type), POINTER :: harmonics
2230 : INTEGER, INTENT(IN) :: nspins
2231 :
2232 : CHARACTER(len=*), PARAMETER :: routineN = 'gaVxcgb_noGC'
2233 :
2234 : INTEGER :: handle, ia, ic, icg, ipgf1, ipgf2, ir, iset1, iset2, iso, iso1, iso2, ispin, l, &
2235 : ld, lmax12, lmax_expansion, lmin12, m1, m2, max_iso_not0, max_iso_not0_local, max_s_harm, &
2236 : maxl, maxso, n1, n2, na, ngau1, ngau2, nngau1, nr, nset, size1
2237 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
2238 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
2239 11965 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2240 11965 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
2241 11965 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gg, gVg_h, gVg_s, matso_h, matso_s, vx
2242 11965 : REAL(dp), DIMENSION(:, :), POINTER :: zet
2243 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
2244 :
2245 11965 : CALL timeset(routineN, handle)
2246 :
2247 11965 : NULLIFY (lmin, lmax, npgf, zet, my_CG)
2248 :
2249 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
2250 : maxso=maxso, maxl=maxl, npgf=npgf, &
2251 11965 : nset=nset, zet=zet)
2252 :
2253 11965 : nr = grid_atom%nr
2254 11965 : na = grid_atom%ng_sphere
2255 11965 : my_CG => harmonics%my_CG
2256 11965 : max_iso_not0 = harmonics%max_iso_not0
2257 11965 : lmax_expansion = indso(1, max_iso_not0)
2258 11965 : max_s_harm = harmonics%max_s_harm
2259 :
2260 83755 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
2261 71790 : ALLOCATE (gVg_h(na, 0:2*maxl), gVg_s(na, 0:2*maxl))
2262 : ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
2263 71790 : matso_s(nsoset(maxl), nsoset(maxl)))
2264 47860 : ALLOCATE (vx(na, nr))
2265 71790 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
2266 :
2267 756515 : g1 = 0.0_dp
2268 756515 : g2 = 0.0_dp
2269 11965 : m1 = 0
2270 40316 : DO iset1 = 1, nset
2271 28351 : n1 = nsoset(lmax(iset1))
2272 28351 : m2 = 0
2273 115582 : DO iset2 = 1, nset
2274 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2275 87231 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2276 87231 : CPASSERT(max_iso_not0_local <= max_iso_not0)
2277 :
2278 87231 : n2 = nsoset(lmax(iset2))
2279 314766 : DO ipgf1 = 1, npgf(iset1)
2280 227535 : ngau1 = n1*(ipgf1 - 1) + m1
2281 227535 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2282 227535 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2283 :
2284 13728285 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2285 1067498 : DO ipgf2 = 1, npgf(iset2)
2286 752732 : ngau2 = n2*(ipgf2 - 1) + m2
2287 :
2288 44880032 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2289 752732 : lmin12 = lmin(iset1) + lmin(iset2)
2290 752732 : lmax12 = lmax(iset1) + lmax(iset2)
2291 :
2292 : ! reduce expansion local densities
2293 980267 : IF (lmin12 <= lmax_expansion) THEN
2294 :
2295 199994374 : gg = 0.0_dp
2296 751787 : IF (lmin12 == 0) THEN
2297 26366556 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2298 : ELSE
2299 18465281 : gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
2300 : END IF
2301 :
2302 : ! limit the expansion of the local densities to a max L
2303 751787 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
2304 :
2305 1067015 : DO l = lmin12 + 1, lmax12
2306 21329215 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
2307 : END DO
2308 :
2309 1658498 : DO ispin = 1, nspins
2310 906711 : ld = lmax12 + 1
2311 57745261 : DO ir = 1, nr
2312 2899672761 : vx(1:na, ir) = vxc_h(1:na, ir, ispin)
2313 : END DO
2314 : CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, vx(1:na, 1:nr), na, &
2315 906711 : gg(1:nr, 0:lmax12), nr, 0.0_dp, gVg_h(1:na, 0:lmax12), na)
2316 57745261 : DO ir = 1, nr
2317 2899672761 : vx(1:na, ir) = vxc_s(1:na, ir, ispin)
2318 : END DO
2319 : CALL dgemm('N', 'N', na, ld, nr, 1.0_dp, vx(1:na, 1:nr), na, &
2320 906711 : gg(1:nr, 0:lmax12), nr, 0.0_dp, gVg_s(1:na, 0:lmax12), na)
2321 :
2322 84907391 : matso_h = 0.0_dp
2323 84907391 : matso_s = 0.0_dp
2324 6961056 : DO iso = 1, max_iso_not0_local
2325 18514785 : DO icg = 1, cg_n_list(iso)
2326 11553729 : iso1 = cg_list(1, icg, iso)
2327 11553729 : iso2 = cg_list(2, icg, iso)
2328 11553729 : l = indso(1, iso1) + indso(1, iso2)
2329 :
2330 11553729 : CPASSERT(l <= lmax_expansion)
2331 595294524 : DO ia = 1, na
2332 : matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2333 : gVg_h(ia, l)* &
2334 : my_CG(iso1, iso2, iso)* &
2335 577686450 : harmonics%slm(ia, iso)
2336 : matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2337 : gVg_s(ia, l)* &
2338 : my_CG(iso1, iso2, iso)* &
2339 589240179 : harmonics%slm(ia, iso)
2340 : END DO
2341 : END DO
2342 : END DO
2343 :
2344 : ! Write in the global matrix
2345 3957549 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2346 2299051 : iso1 = nsoset(lmin(iset1) - 1) + 1
2347 2299051 : iso2 = ngau2 + ic
2348 : CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2349 2299051 : int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2350 : CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2351 3205762 : int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
2352 : END DO
2353 :
2354 : END DO ! ispin
2355 :
2356 : END IF ! lmax_expansion
2357 :
2358 : END DO ! ipfg2
2359 : END DO ! ipfg1
2360 202813 : m2 = m2 + maxso
2361 : END DO ! iset2
2362 40316 : m1 = m1 + maxso
2363 : END DO ! iset1
2364 :
2365 11965 : DEALLOCATE (g1, g2, gg, matso_h, matso_s, gVg_s, gVg_h, vx)
2366 :
2367 11965 : DEALLOCATE (cg_list, cg_n_list)
2368 :
2369 11965 : CALL timestop(handle)
2370 :
2371 11965 : END SUBROUTINE gaVxcgb_noGC
2372 :
2373 : ! **************************************************************************************************
2374 : !> \brief ...
2375 : !> \param vxc_h ...
2376 : !> \param vxc_s ...
2377 : !> \param vxg_h ...
2378 : !> \param vxg_s ...
2379 : !> \param int_hh ...
2380 : !> \param int_ss ...
2381 : !> \param grid_atom ...
2382 : !> \param basis_1c ...
2383 : !> \param harmonics ...
2384 : !> \param nspins ...
2385 : ! **************************************************************************************************
2386 18211 : SUBROUTINE gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
2387 : grid_atom, basis_1c, harmonics, nspins)
2388 :
2389 : REAL(dp), DIMENSION(:, :, :), POINTER :: vxc_h, vxc_s
2390 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: vxg_h, vxg_s
2391 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
2392 : TYPE(grid_atom_type), POINTER :: grid_atom
2393 : TYPE(gto_basis_set_type), POINTER :: basis_1c
2394 : TYPE(harmonics_atom_type), POINTER :: harmonics
2395 : INTEGER, INTENT(IN) :: nspins
2396 :
2397 : CHARACTER(len=*), PARAMETER :: routineN = 'gaVxcgb_GC'
2398 :
2399 : INTEGER :: dmax_iso_not0_local, handle, ia, ic, icg, ipgf1, ipgf2, ir, iset1, iset2, iso, &
2400 : iso1, iso2, ispin, l, lmax12, lmax_expansion, lmin12, m1, m2, max_iso_not0, &
2401 : max_iso_not0_local, max_s_harm, maxl, maxso, n1, n2, na, ngau1, ngau2, nngau1, nr, nset, &
2402 : size1
2403 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list, dcg_n_list
2404 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list, dcg_list
2405 18211 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2406 : REAL(dp) :: urad
2407 18211 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
2408 18211 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gVXCg_h, gVXCg_s, matso_h, &
2409 18211 : matso_s
2410 18211 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: gVXGg_h, gVXGg_s
2411 18211 : REAL(dp), DIMENSION(:, :), POINTER :: zet
2412 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
2413 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz
2414 :
2415 18211 : CALL timeset(routineN, handle)
2416 :
2417 18211 : NULLIFY (lmin, lmax, npgf, zet, my_CG, my_CG_dxyz)
2418 :
2419 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
2420 : maxso=maxso, maxl=maxl, npgf=npgf, &
2421 18211 : nset=nset, zet=zet)
2422 :
2423 18211 : nr = grid_atom%nr
2424 18211 : na = grid_atom%ng_sphere
2425 18211 : my_CG => harmonics%my_CG
2426 18211 : my_CG_dxyz => harmonics%my_CG_dxyz
2427 18211 : max_iso_not0 = harmonics%max_iso_not0
2428 18211 : lmax_expansion = indso(1, max_iso_not0)
2429 18211 : max_s_harm = harmonics%max_s_harm
2430 :
2431 163899 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl))
2432 109266 : ALLOCATE (gVXCg_h(na, 0:2*maxl), gVXCg_s(na, 0:2*maxl))
2433 109266 : ALLOCATE (gVXGg_h(3, na, 0:2*maxl), gVXGg_s(3, na, 0:2*maxl))
2434 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm), &
2435 163899 : dcg_list(2, nsoset(maxl)**2, max_s_harm), dcg_n_list(max_s_harm))
2436 :
2437 : ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
2438 109266 : matso_s(nsoset(maxl), nsoset(maxl)))
2439 :
2440 39040 : DO ispin = 1, nspins
2441 :
2442 1087359 : g1 = 0.0_dp
2443 1087359 : g2 = 0.0_dp
2444 20829 : m1 = 0
2445 91131 : DO iset1 = 1, nset
2446 52091 : n1 = nsoset(lmax(iset1))
2447 52091 : m2 = 0
2448 231156 : DO iset2 = 1, nset
2449 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2450 179065 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2451 179065 : CPASSERT(max_iso_not0_local <= max_iso_not0)
2452 : CALL get_none0_cg_list(my_CG_dxyz, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2453 179065 : max_s_harm, lmax_expansion, dcg_list, dcg_n_list, dmax_iso_not0_local)
2454 :
2455 179065 : n2 = nsoset(lmax(iset2))
2456 558628 : DO ipgf1 = 1, npgf(iset1)
2457 379563 : ngau1 = n1*(ipgf1 - 1) + m1
2458 379563 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2459 379563 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2460 :
2461 20078793 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2462 1474193 : DO ipgf2 = 1, npgf(iset2)
2463 915565 : ngau2 = n2*(ipgf2 - 1) + m2
2464 :
2465 48720495 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2466 915565 : lmin12 = lmin(iset1) + lmin(iset2)
2467 915565 : lmax12 = lmax(iset1) + lmax(iset2)
2468 :
2469 : !test reduce expansion local densities
2470 915565 : IF (lmin12 <= lmax_expansion) THEN
2471 :
2472 205545210 : gg = 0.0_dp
2473 205545210 : dgg = 0.0_dp
2474 :
2475 914965 : IF (lmin12 == 0) THEN
2476 32951634 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2477 : ELSE
2478 15738261 : gg(1:nr, lmin12) = grid_atom%rad2l(1:nr, lmin12)*g1(1:nr)*g2(1:nr)
2479 : END IF
2480 :
2481 : !test reduce expansion local densities
2482 914965 : IF (lmax12 > lmax_expansion) lmax12 = lmax_expansion
2483 :
2484 1434285 : DO l = lmin12 + 1, lmax12
2485 27944120 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
2486 : dgg(1:nr, l - 1) = dgg(1:nr, l - 1) - 2.0_dp*(zet(ipgf1, iset1) + &
2487 28859085 : zet(ipgf2, iset2))*gg(1:nr, l)
2488 : END DO
2489 : dgg(1:nr, lmax12) = dgg(1:nr, lmax12) - 2.0_dp*(zet(ipgf1, iset1) + &
2490 : zet(ipgf2, iset2))*grid_atom%rad(1:nr)* &
2491 48689895 : gg(1:nr, lmax12)
2492 :
2493 195781486 : gVXCg_h = 0.0_dp
2494 195781486 : gVXCg_s = 0.0_dp
2495 768909544 : gVXGg_h = 0.0_dp
2496 768909544 : gVXGg_s = 0.0_dp
2497 :
2498 : ! Cross Term
2499 2349250 : DO l = lmin12, lmax12
2500 74014732 : DO ia = 1, na
2501 3836628667 : DO ir = 1, nr
2502 : gVXCg_h(ia, l) = gVXCg_h(ia, l) + &
2503 : gg(ir, l)*vxc_h(ia, ir, ispin) + &
2504 : dgg(ir, l)* &
2505 : (vxg_h(1, ia, ir, ispin)*harmonics%a(1, ia) + &
2506 : vxg_h(2, ia, ir, ispin)*harmonics%a(2, ia) + &
2507 3763528900 : vxg_h(3, ia, ir, ispin)*harmonics%a(3, ia))
2508 :
2509 : gVXCg_s(ia, l) = gVXCg_s(ia, l) + &
2510 : gg(ir, l)*vxc_s(ia, ir, ispin) + &
2511 : dgg(ir, l)* &
2512 : (vxg_s(1, ia, ir, ispin)*harmonics%a(1, ia) + &
2513 : vxg_s(2, ia, ir, ispin)*harmonics%a(2, ia) + &
2514 3763528900 : vxg_s(3, ia, ir, ispin)*harmonics%a(3, ia))
2515 :
2516 3763528900 : urad = grid_atom%oorad2l(ir, 1)
2517 :
2518 : gVXGg_h(1, ia, l) = gVXGg_h(1, ia, l) + &
2519 : vxg_h(1, ia, ir, ispin)* &
2520 3763528900 : gg(ir, l)*urad
2521 :
2522 : gVXGg_h(2, ia, l) = gVXGg_h(2, ia, l) + &
2523 : vxg_h(2, ia, ir, ispin)* &
2524 3763528900 : gg(ir, l)*urad
2525 :
2526 : gVXGg_h(3, ia, l) = gVXGg_h(3, ia, l) + &
2527 : vxg_h(3, ia, ir, ispin)* &
2528 3763528900 : gg(ir, l)*urad
2529 :
2530 : gVXGg_s(1, ia, l) = gVXGg_s(1, ia, l) + &
2531 : vxg_s(1, ia, ir, ispin)* &
2532 3763528900 : gg(ir, l)*urad
2533 :
2534 : gVXGg_s(2, ia, l) = gVXGg_s(2, ia, l) + &
2535 : vxg_s(2, ia, ir, ispin)* &
2536 3763528900 : gg(ir, l)*urad
2537 :
2538 : gVXGg_s(3, ia, l) = gVXGg_s(3, ia, l) + &
2539 : vxg_s(3, ia, ir, ispin)* &
2540 3835194382 : gg(ir, l)*urad
2541 :
2542 : END DO ! ir
2543 : END DO ! ia
2544 : END DO ! l
2545 :
2546 72488613 : matso_h = 0.0_dp
2547 72488613 : matso_s = 0.0_dp
2548 6079146 : DO iso = 1, max_iso_not0_local
2549 16577599 : DO icg = 1, cg_n_list(iso)
2550 10498453 : iso1 = cg_list(1, icg, iso)
2551 10498453 : iso2 = cg_list(2, icg, iso)
2552 :
2553 10498453 : l = indso(1, iso1) + indso(1, iso2)
2554 :
2555 : !test reduce expansion local densities
2556 10498453 : CPASSERT(l <= lmax_expansion)
2557 540336944 : DO ia = 1, na
2558 : matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2559 : gVXCg_h(ia, l)* &
2560 : harmonics%slm(ia, iso)* &
2561 524674310 : my_CG(iso1, iso2, iso)
2562 : matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2563 : gVXCg_s(ia, l)* &
2564 : harmonics%slm(ia, iso)* &
2565 535172763 : my_CG(iso1, iso2, iso)
2566 : END DO ! ia
2567 :
2568 : !test reduce expansion local densities
2569 :
2570 : END DO
2571 :
2572 : END DO ! iso
2573 :
2574 3229265 : DO iso = 1, dmax_iso_not0_local
2575 19387957 : DO icg = 1, dcg_n_list(iso)
2576 16158692 : iso1 = dcg_list(1, icg, iso)
2577 16158692 : iso2 = dcg_list(2, icg, iso)
2578 :
2579 16158692 : l = indso(1, iso1) + indso(1, iso2)
2580 : !test reduce expansion local densities
2581 16158692 : CPASSERT(l <= lmax_expansion)
2582 825668260 : DO ia = 1, na
2583 : matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2584 : (gVXGg_h(1, ia, l)*my_CG_dxyz(1, iso1, iso2, iso) + &
2585 : gVXGg_h(2, ia, l)*my_CG_dxyz(2, iso1, iso2, iso) + &
2586 : gVXGg_h(3, ia, l)*my_CG_dxyz(3, iso1, iso2, iso))* &
2587 807195268 : harmonics%slm(ia, iso)
2588 :
2589 : matso_s(iso1, iso2) = matso_s(iso1, iso2) + &
2590 : (gVXGg_s(1, ia, l)*my_CG_dxyz(1, iso1, iso2, iso) + &
2591 : gVXGg_s(2, ia, l)*my_CG_dxyz(2, iso1, iso2, iso) + &
2592 : gVXGg_s(3, ia, l)*my_CG_dxyz(3, iso1, iso2, iso))* &
2593 823353960 : harmonics%slm(ia, iso)
2594 :
2595 : END DO ! ia
2596 :
2597 : !test reduce expansion local densities
2598 :
2599 : END DO ! icg
2600 : END DO ! iso
2601 : !test reduce expansion local densities
2602 : END IF ! lmax_expansion
2603 :
2604 : ! Write in the global matrix
2605 3520997 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2606 2225869 : iso1 = nsoset(lmin(iset1) - 1) + 1
2607 2225869 : iso2 = ngau2 + ic
2608 : CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2609 2225869 : int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2610 : CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2611 3141434 : int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
2612 : END DO
2613 :
2614 : END DO ! ipfg2
2615 : END DO ! ipfg1
2616 589286 : m2 = m2 + maxso
2617 : END DO ! iset2
2618 72920 : m1 = m1 + maxso
2619 : END DO ! iset1
2620 : END DO ! ispin
2621 :
2622 18211 : DEALLOCATE (g1, g2, gg, dgg, matso_h, matso_s, gVXCg_h, gVXCg_s, gVXGg_h, gVXGg_s)
2623 18211 : DEALLOCATE (cg_list, cg_n_list, dcg_list, dcg_n_list)
2624 :
2625 18211 : CALL timestop(handle)
2626 :
2627 18211 : END SUBROUTINE gaVxcgb_GC
2628 :
2629 : ! **************************************************************************************************
2630 : !> \brief Integrates 0.5 * grad_ga .dot. (V_tau * grad_gb) on the atomic grid for meta-GGA
2631 : !> \param vtau_h the har tau potential
2632 : !> \param vtau_s the soft tau potential
2633 : !> \param int_hh ...
2634 : !> \param int_ss ...
2635 : !> \param grid_atom ...
2636 : !> \param basis_1c ...
2637 : !> \param harmonics ...
2638 : !> \param nspins ...
2639 : !> \note This is a rewrite to correct meta-GGA GAPW bug. This is more brute force than the original
2640 : !> but makes sure that no corner is cut in terms of accuracy (A. Bussy)
2641 : ! **************************************************************************************************
2642 356 : SUBROUTINE dgaVtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
2643 : grid_atom, basis_1c, harmonics, nspins)
2644 :
2645 : REAL(dp), DIMENSION(:, :, :), POINTER :: vtau_h, vtau_s
2646 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_hh, int_ss
2647 : TYPE(grid_atom_type), POINTER :: grid_atom
2648 : TYPE(gto_basis_set_type), POINTER :: basis_1c
2649 : TYPE(harmonics_atom_type), POINTER :: harmonics
2650 : INTEGER, INTENT(IN) :: nspins
2651 :
2652 : CHARACTER(len=*), PARAMETER :: routineN = 'dgaVtaudgb'
2653 :
2654 : INTEGER :: dir, handle, ipgf, iset, iso, ispin, &
2655 : jpgf, jset, jso, l, maxso, na, nr, &
2656 : nset, starti, startj
2657 356 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2658 356 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fga, fgr, r1, r2, work
2659 356 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: a1, a2, intso_h, intso_s
2660 356 : REAL(dp), DIMENSION(:, :), POINTER :: slm, zet
2661 356 : REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
2662 :
2663 356 : CALL timeset(routineN, handle)
2664 :
2665 356 : NULLIFY (zet, slm, dslm_dxyz, lmax, lmin, npgf)
2666 :
2667 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
2668 356 : maxso=maxso, npgf=npgf, nset=nset, zet=zet)
2669 :
2670 356 : na = grid_atom%ng_sphere
2671 356 : nr = grid_atom%nr
2672 :
2673 356 : slm => harmonics%slm
2674 356 : dslm_dxyz => harmonics%dslm_dxyz
2675 :
2676 : ! Separate the functions into purely r and purely angular parts and precompute them all
2677 : ! Not memory intensive since 1D arrays
2678 2492 : ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3))
2679 2136 : ALLOCATE (r1(nr, nset*maxso), r2(nr, nset*maxso))
2680 4521212 : a1 = 0.0_dp; a2 = 0.0_dp
2681 1521156 : r1 = 0.0_dp; r2 = 0.0_dp
2682 :
2683 1168 : DO iset = 1, nset
2684 3632 : DO ipgf = 1, npgf(iset)
2685 2464 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
2686 10307 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
2687 7031 : l = indso(1, iso)
2688 :
2689 : ! The x derivative of the spherical orbital, divided in angular and radial parts
2690 : ! Two of each are needed because d/dx(r^l Y_lm) * exp(-al*r^2) + r^l Y_lm * d/dx(exp-al*r^2)
2691 :
2692 : ! the purely radial part of d/dx(r^l Y_lm) * exp(-al*r^2) (same for y,z)
2693 376781 : r1(1:nr, starti + iso) = grid_atom%rad(1:nr)**(l - 1)*EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
2694 :
2695 : ! the purely radial part of r^l Y_lm * d/dx(exp-al*r^2) (same for y,z)
2696 : r2(1:nr, starti + iso) = -2.0_dp*zet(ipgf, iset)*grid_atom%rad(1:nr)**(l + 1) &
2697 376781 : *EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
2698 :
2699 30588 : DO dir = 1, 3
2700 : ! the purely angular part of d/dx(r^l Y_lm) * exp(-al*r^2)
2701 1115055 : a1(1:na, starti + iso, dir) = dslm_dxyz(dir, 1:na, iso)
2702 :
2703 : ! the purely angular part of r^l Y_lm * d/dx(exp-al*r^2)
2704 1122086 : a2(1:na, starti + iso, dir) = harmonics%a(dir, 1:na)*slm(1:na, iso)
2705 : END DO
2706 :
2707 : END DO !iso
2708 : END DO !ipgf
2709 : END DO !iset
2710 :
2711 : !Do the integration in terms of matrix-matrix multiplications
2712 1780 : ALLOCATE (intso_h(nset*maxso, nset*maxso, nspins))
2713 1424 : ALLOCATE (intso_s(nset*maxso, nset*maxso, nspins))
2714 5310320 : intso_h = 0.0_dp; intso_s = 0.0_dp
2715 :
2716 1068 : ALLOCATE (fga(na, 1))
2717 1068 : ALLOCATE (fgr(nr, 1))
2718 712 : ALLOCATE (work(na, 1))
2719 57996 : fga = 0.0_dp; fgr = 0.0_dp; work = 0.0_dp
2720 :
2721 1168 : DO iset = 1, nset
2722 3576 : DO jset = 1, nset
2723 10518 : DO ipgf = 1, npgf(iset)
2724 7298 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
2725 36446 : DO jpgf = 1, npgf(jset)
2726 26740 : startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
2727 :
2728 110858 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
2729 376637 : DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
2730 :
2731 : !Two component per function => 4 terms in total
2732 :
2733 : ! take 0.5*r1*a1(dir) * V * r1*a1(dir)
2734 14245427 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r1(1:nr, startj + jso)
2735 1092308 : DO dir = 1, 3
2736 42468741 : fga(1:na, 1) = a1(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
2737 :
2738 1911539 : DO ispin = 1, nspins
2739 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2740 819231 : nr, 0.0_dp, work, na)
2741 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2742 819231 : intso_h(starti + iso:, startj + jso, ispin), 1)
2743 :
2744 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_s(:, :, ispin), na, fgr, &
2745 819231 : nr, 0.0_dp, work, na)
2746 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2747 1638462 : intso_s(starti + iso:, startj + jso, ispin), 1)
2748 : END DO
2749 : END DO !dir
2750 :
2751 : ! add 0.5*r1*a1(dir) * V * r2*a2(dir)
2752 14245427 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r2(1:nr, startj + jso)
2753 1092308 : DO dir = 1, 3
2754 42468741 : fga(1:na, 1) = a1(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
2755 :
2756 1911539 : DO ispin = 1, nspins
2757 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2758 819231 : nr, 0.0_dp, work, na)
2759 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2760 819231 : intso_h(starti + iso:, startj + jso, ispin), 1)
2761 :
2762 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_s(:, :, ispin), na, fgr, &
2763 819231 : nr, 0.0_dp, work, na)
2764 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2765 1638462 : intso_s(starti + iso:, startj + jso, ispin), 1)
2766 : END DO
2767 : END DO !dir
2768 :
2769 : ! add 0.5*r2*a2(dir) * V * r1*a1(dir)
2770 14245427 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r1(1:nr, startj + jso)
2771 1092308 : DO dir = 1, 3
2772 42468741 : fga(1:na, 1) = a2(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
2773 :
2774 1911539 : DO ispin = 1, nspins
2775 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2776 819231 : nr, 0.0_dp, work, na)
2777 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2778 819231 : intso_h(starti + iso:, startj + jso, ispin), 1)
2779 :
2780 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_s(:, :, ispin), na, fgr, &
2781 819231 : nr, 0.0_dp, work, na)
2782 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2783 1638462 : intso_s(starti + iso:, startj + jso, ispin), 1)
2784 : END DO
2785 : END DO !dir
2786 :
2787 : ! add the last term: 0.5*r2*a2(dir) * V * r2*a2(dir)
2788 14245427 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r2(1:nr, startj + jso)
2789 1169128 : DO dir = 1, 3
2790 42468741 : fga(1:na, 1) = a2(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
2791 :
2792 1911539 : DO ispin = 1, nspins
2793 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2794 819231 : nr, 0.0_dp, work, na)
2795 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2796 819231 : intso_h(starti + iso:, startj + jso, ispin), 1)
2797 :
2798 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_s(:, :, ispin), na, fgr, &
2799 819231 : nr, 0.0_dp, work, na)
2800 : CALL dgemm('T', 'N', 1, 1, na, 1.0_dp, work, na, fga, na, 1.0_dp, &
2801 1638462 : intso_s(starti + iso:, startj + jso, ispin), 1)
2802 : END DO
2803 : END DO !dir
2804 :
2805 : END DO !jso
2806 : END DO !iso
2807 :
2808 : END DO !jpgf
2809 : END DO !ipgf
2810 : END DO !jset
2811 : END DO !iset
2812 :
2813 : ! Put the integrals in the rho_atom data structure
2814 712 : DO ispin = 1, nspins
2815 2654982 : int_hh(ispin)%r_coef(:, :) = int_hh(ispin)%r_coef(:, :) + intso_h(:, :, ispin)
2816 2655338 : int_ss(ispin)%r_coef(:, :) = int_ss(ispin)%r_coef(:, :) + intso_s(:, :, ispin)
2817 : END DO
2818 :
2819 356 : CALL timestop(handle)
2820 :
2821 712 : END SUBROUTINE dgaVtaudgb
2822 :
2823 : END MODULE qs_vxc_atom
|
