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