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