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