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