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 57162 : 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 19054 : 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 38108 : REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight
112 38108 : REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
113 19054 : vtau_s, vxc_h, vxc_s
114 38108 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
115 19054 : 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 19054 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
122 19054 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
123 19054 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
124 19054 : 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 19054 : CALL timeset(routineN, handle)
134 :
135 19054 : NULLIFY (atom_list)
136 19054 : NULLIFY (my_kind_set)
137 19054 : NULLIFY (atomic_kind_set)
138 19054 : NULLIFY (grid_atom)
139 19054 : NULLIFY (harmonics)
140 19054 : NULLIFY (input)
141 19054 : NULLIFY (para_env)
142 19054 : NULLIFY (rho_atom)
143 19054 : NULLIFY (my_rho_atom_set)
144 19054 : NULLIFY (rho_nlcc)
145 :
146 19054 : epr_xc = .FALSE.
147 19054 : IF (PRESENT(gradient_atom_set)) THEN
148 10 : epr_xc = .TRUE.
149 : END IF
150 :
151 19054 : IF (PRESENT(adiabatic_rescale_factor)) THEN
152 44 : my_adiabatic_rescale_factor = adiabatic_rescale_factor
153 : ELSE
154 19010 : 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 19054 : rho_atom_set=my_rho_atom_set)
164 :
165 19054 : IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
166 19054 : IF (PRESENT(rho_atom_set_external)) my_rho_atom_set => rho_atom_set_external
167 :
168 19054 : nlcc = has_nlcc(my_kind_set)
169 :
170 19054 : 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 19044 : my_xc_section => section_vals_get_subs_vals(input, "DFT%XC")
175 : END IF
176 :
177 19054 : IF (PRESENT(xc_section_external)) my_xc_section => xc_section_external
178 :
179 19054 : 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 19054 : i_val=myfun)
182 :
183 19054 : 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 16614 : r_val=density_cut)
190 : CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
191 16614 : r_val=gradient_cut)
192 : CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
193 16614 : r_val=tau_cut)
194 :
195 16614 : lsd = dft_control%lsd
196 16614 : nspins = dft_control%nspins
197 : needs = xc_functionals_get_needs(xc_fun_section, &
198 : lsd=lsd, &
199 16614 : calc_potential=.TRUE.)
200 :
201 : ! whatever the xc, if epr_xc, drho_spin is needed
202 16614 : IF (epr_xc) needs%drho_spin = .TRUE.
203 :
204 16614 : gradient_f = (needs%drho .OR. needs%drho_spin)
205 16614 : tau_f = (needs%tau .OR. needs%tau_spin)
206 :
207 : ! Initialize energy contribution from the one center XC terms to zero
208 16614 : exc1 = 0.0_dp
209 :
210 : ! Nullify some pointers for work-arrays
211 16614 : NULLIFY (rho_h, drho_h, rho_s, drho_s, weight)
212 16614 : NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
213 16614 : NULLIFY (tau_h, tau_s)
214 16614 : NULLIFY (vtau_h, vtau_s)
215 :
216 : ! Here starts the loop over all the atoms
217 :
218 49960 : DO ikind = 1, SIZE(atomic_kind_set)
219 33346 : 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 33346 : harmonics=harmonics, grid_atom=grid_atom)
222 33346 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
223 :
224 33346 : IF (.NOT. paw_atom) CYCLE
225 :
226 30538 : nr = grid_atom%nr
227 30538 : 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 305380 : bounds(1:2, 1:3) = 1
235 30538 : bounds(2, 1) = na
236 30538 : bounds(2, 2) = nr
237 :
238 : ! create a place where to put the derivatives
239 30538 : 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 30538 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
243 : CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
244 30538 : drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
245 :
246 : ! allocate the required 3d arrays where to store rho and drho
247 30538 : CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
248 30538 : CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
249 :
250 30538 : CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
251 30538 : CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
252 30538 : weight => grid_atom%weight
253 30538 : CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
254 30538 : CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
255 : !
256 30538 : IF (gradient_f) THEN
257 18308 : CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
258 18308 : CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
259 18308 : CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
260 18308 : CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
261 : END IF
262 :
263 30538 : IF (tau_f) THEN
264 444 : CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
265 444 : CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
266 444 : CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
267 444 : 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 30538 : donlcc = .FALSE.
272 30538 : 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 30538 : num_pe = para_env%num_pe
281 30538 : bo = get_limit(natom, para_env%num_pe, para_env%mepos)
282 :
283 54940 : DO iat = bo(1), bo(2)
284 24402 : iatom = atom_list(iat)
285 :
286 24402 : my_rho_atom_set(iatom)%exc_h = 0.0_dp
287 24402 : my_rho_atom_set(iatom)%exc_s = 0.0_dp
288 :
289 24402 : rho_atom => my_rho_atom_set(iatom)
290 86040274 : rho_h = 0.0_dp
291 86040274 : rho_s = 0.0_dp
292 24402 : IF (gradient_f) THEN
293 14369 : 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 14369 : rho_rad_s_d=r_s_d)
298 222077009 : drho_h = 0.0_dp
299 222077009 : drho_s = 0.0_dp
300 : ELSE
301 10033 : NULLIFY (r_h, r_s)
302 10033 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
303 10033 : rho_d = 0.0_dp
304 : END IF
305 24402 : IF (tau_f) THEN
306 : !compute tau on the grid all at once
307 308 : CALL calc_tau_atom(tau_h, tau_s, rho_atom, my_kind_set(ikind), nspins)
308 : ELSE
309 24094 : tau_d = 0.0_dp
310 : END IF
311 :
312 1415042 : 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 1390640 : r_h_d, r_s_d, drho_h, drho_s)
316 1415042 : 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 1415042 : DO ir = 1, nr
322 1415042 : IF (tau_f) THEN
323 16700 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
324 16700 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
325 1373940 : ELSE IF (gradient_f) THEN
326 719190 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
327 719190 : CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
328 : ELSE
329 654750 : CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
330 654750 : 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 24402 : 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 24402 : epr_xc=epr_xc, adiabatic_rescale_factor=my_adiabatic_rescale_factor)
341 24402 : rho_atom%exc_h = rho_atom%exc_h + exc_h
342 :
343 : !-------------------!
344 : ! soft atom density !
345 : !-------------------!
346 24402 : 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 24402 : epr_xc=epr_xc, adiabatic_rescale_factor=my_adiabatic_rescale_factor)
350 24402 : rho_atom%exc_s = rho_atom%exc_s + exc_s
351 :
352 24402 : 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 24402 : 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 24402 : IF (.NOT. energy_only) THEN
378 22980 : NULLIFY (int_hh, int_ss)
379 22980 : CALL get_rho_atom(rho_atom=rho_atom, ga_Vlocal_gb_h=int_hh, ga_Vlocal_gb_s=int_ss)
380 22980 : IF (gradient_f) THEN
381 : CALL gaVxcgb_GC(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
382 13129 : grid_atom, basis_1c, harmonics, nspins)
383 : ELSE
384 : CALL gaVxcgb_noGC(vxc_h, vxc_s, int_hh, int_ss, &
385 9851 : grid_atom, basis_1c, harmonics, nspins)
386 : END IF
387 22980 : IF (tau_f) THEN
388 : CALL dgaVtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
389 308 : grid_atom, basis_1c, harmonics, nspins)
390 : END IF
391 : END IF ! energy_only
392 54940 : 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 30538 : CALL xc_dset_release(deriv_set)
397 30538 : CALL xc_rho_set_release(rho_set_h)
398 111036 : CALL xc_rho_set_release(rho_set_s)
399 :
400 : END DO ! ikind
401 :
402 16614 : CALL para_env%sum(exc1)
403 :
404 16614 : IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
405 16614 : IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
406 16614 : IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
407 16614 : IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
408 :
409 16614 : IF (gradient_f) THEN
410 10098 : IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
411 10098 : IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
412 10098 : IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
413 10098 : IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
414 : END IF
415 :
416 16614 : IF (tau_f) THEN
417 244 : IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
418 244 : IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
419 244 : IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
420 244 : IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
421 : END IF
422 :
423 : END IF !xc_none
424 :
425 19054 : CALL timestop(handle)
426 :
427 819322 : 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 1819440 : 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 1819440 : CPASSERT(ASSOCIATED(r_h))
1533 1819440 : CPASSERT(ASSOCIATED(r_s))
1534 1819440 : CPASSERT(ASSOCIATED(rho_h))
1535 1819440 : CPASSERT(ASSOCIATED(rho_s))
1536 1819440 : IF (grad_func) THEN
1537 1053590 : CPASSERT(ASSOCIATED(dr_h))
1538 1053590 : CPASSERT(ASSOCIATED(dr_s))
1539 1053590 : CPASSERT(ASSOCIATED(r_h_d))
1540 1053590 : CPASSERT(ASSOCIATED(r_s_d))
1541 1053590 : CPASSERT(ASSOCIATED(drho_h))
1542 1053590 : CPASSERT(ASSOCIATED(drho_s))
1543 : END IF
1544 :
1545 1819440 : na = grid_atom%ng_sphere
1546 1819440 : rad = grid_atom%rad(ir)
1547 1819440 : urad = grid_atom%oorad2l(ir, 1)
1548 3939120 : DO ispin = 1, nspins
1549 31337040 : DO iso = 1, harmonics%max_iso_not0
1550 1400955120 : DO ia = 1, na
1551 : rho_h(ia, ir, ispin) = rho_h(ia, ir, ispin) + &
1552 1371437520 : r_h(ispin)%r_coef(ir, iso)*harmonics%slm(ia, iso)
1553 : rho_s(ia, ir, ispin) = rho_s(ia, ir, ispin) + &
1554 1398835440 : 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 1819440 : IF (grad_func) THEN
1560 2260770 : DO ispin = 1, nspins
1561 18552990 : DO iso = 1, harmonics%max_iso_not0
1562 833651920 : 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 816152520 : 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 816152520 : 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 816152520 : 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 816152520 : 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 816152520 : 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 816152520 : 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 816152520 : 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 832444740 : 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 1819440 : 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 308 : 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 308 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, o2nindex
1643 : REAL(dp) :: cpc_h, cpc_s
1644 308 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fga, fgr, r1, r2
1645 308 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: a1, a2
1646 308 : REAL(dp), DIMENSION(:, :), POINTER :: slm, zet
1647 308 : 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 308 : NULLIFY (harmonics, grid_atom, basis_1c, lmax, lmin, npgf, zet, slm, dslm_dxyz, o2nindex)
1653 :
1654 308 : 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 308 : CALL get_qs_kind(qs_kind, grid_atom=grid_atom, harmonics=harmonics)
1659 308 : CALL get_qs_kind(qs_kind, basis_set=basis_1c, basis_type="GAPW_1C")
1660 :
1661 308 : nr = grid_atom%nr
1662 308 : na = grid_atom%ng_sphere
1663 :
1664 308 : slm => harmonics%slm
1665 308 : dslm_dxyz => harmonics%dslm_dxyz
1666 :
1667 308 : CALL get_paw_basis_info(basis_1c, o2nindex=o2nindex)
1668 :
1669 : !zeroing tau, assuming it is already allocated
1670 945916 : tau_h = 0.0_dp
1671 945916 : tau_s = 0.0_dp
1672 :
1673 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, npgf=npgf, &
1674 308 : nset=nset, zet=zet, maxso=maxso)
1675 :
1676 : !Separate the functions into purely r and purely angular parts, precompute them all
1677 2156 : ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3))
1678 1848 : ALLOCATE (r1(nr, nset*maxso), r2(nr, nset*maxso))
1679 3859916 : a1 = 0.0_dp; a2 = 0.0_dp
1680 1300788 : r1 = 0.0_dp; r2 = 0.0_dp
1681 :
1682 1008 : DO iset = 1, nset
1683 3200 : DO ipgf = 1, npgf(iset)
1684 2192 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1685 8211 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
1686 5319 : 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 289469 : 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 289469 : *EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
1697 :
1698 23468 : DO dir = 1, 3
1699 : ! the purely angular part of d/dx(r^l Y_lm) * exp(-al*r^2)
1700 853119 : 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 858438 : 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 924 : ALLOCATE (fga(na, 1))
1712 924 : ALLOCATE (fgr(nr, 1))
1713 33960 : fga = 0.0_dp; fgr = 0.0_dp
1714 :
1715 1008 : DO iset = 1, nset
1716 2984 : DO jset = 1, nset
1717 9254 : DO ipgf = 1, npgf(iset)
1718 6578 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
1719 33710 : DO jpgf = 1, npgf(jset)
1720 25156 : startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
1721 :
1722 99450 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
1723 296525 : DO jso = nsoset(lmin(jset) - 1) + 1, nsoset(lmax(jset))
1724 :
1725 203653 : ip = o2nindex(starti + iso)
1726 203653 : jp = o2nindex(startj + jso)
1727 :
1728 : !Two component per function => 4 terms in total
1729 :
1730 : ! take r1*a1(dir) * r1*a1(dir)
1731 10704803 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r1(1:nr, startj + jso)
1732 814612 : DO dir = 1, 3
1733 31846869 : fga(1:na, 1) = a1(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
1734 :
1735 1425571 : DO ispin = 1, nspins
1736 : !get the projectors
1737 610959 : cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
1738 610959 : cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
1739 :
1740 : !compute contribution to tau
1741 32725368 : DO ir = 1, nr
1742 1676082909 : DO ia = 1, na
1743 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
1744 1643968500 : fgr(ir, 1)*fga(ia, 1)
1745 :
1746 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
1747 1675471950 : 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 10704803 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r2(1:nr, startj + jso)
1756 814612 : DO dir = 1, 3
1757 31846869 : fga(1:na, 1) = a1(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
1758 :
1759 1425571 : DO ispin = 1, nspins
1760 : !get the projectors
1761 610959 : cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
1762 610959 : cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
1763 :
1764 : !compute contribution to tau
1765 32725368 : DO ir = 1, nr
1766 1676082909 : DO ia = 1, na
1767 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
1768 1643968500 : fgr(ir, 1)*fga(ia, 1)
1769 :
1770 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
1771 1675471950 : 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 10704803 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r1(1:nr, startj + jso)
1780 814612 : DO dir = 1, 3
1781 31846869 : fga(1:na, 1) = a2(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
1782 :
1783 1425571 : DO ispin = 1, nspins
1784 : !get the projectors
1785 610959 : cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
1786 610959 : cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
1787 :
1788 : !compute contribution to tau
1789 32725368 : DO ir = 1, nr
1790 1676082909 : DO ia = 1, na
1791 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
1792 1643968500 : fgr(ir, 1)*fga(ia, 1)
1793 :
1794 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
1795 1675471950 : 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 10704803 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r2(1:nr, startj + jso)
1804 882328 : DO dir = 1, 3
1805 31846869 : fga(1:na, 1) = a2(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
1806 :
1807 1425571 : DO ispin = 1, nspins
1808 : !get the projectors
1809 610959 : cpc_h = rho_atom%cpc_h(ispin)%r_coef(ip, jp)
1810 610959 : cpc_s = rho_atom%cpc_s(ispin)%r_coef(ip, jp)
1811 :
1812 : !compute contribution to tau
1813 32725368 : DO ir = 1, nr
1814 1676082909 : DO ia = 1, na
1815 : tau_h(ia, ir, ispin) = tau_h(ia, ir, ispin) + 0.5_dp*cpc_h* &
1816 1643968500 : fgr(ir, 1)*fga(ia, 1)
1817 :
1818 : tau_s(ia, ir, ispin) = tau_s(ia, ir, ispin) + 0.5_dp*cpc_s* &
1819 1675471950 : 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 308 : DEALLOCATE (o2nindex)
1835 :
1836 308 : CALL timestop(handle)
1837 :
1838 924 : 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 11130 : 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 11130 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
1950 11130 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
1951 11130 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: gg, gVg_h, gVg_s, matso_h, matso_s, vx
1952 11130 : REAL(dp), DIMENSION(:, :), POINTER :: zet
1953 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
1954 :
1955 11130 : CALL timeset(routineN, handle)
1956 :
1957 11130 : 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 11130 : nset=nset, zet=zet)
1962 :
1963 11130 : nr = grid_atom%nr
1964 11130 : na = grid_atom%ng_sphere
1965 11130 : my_CG => harmonics%my_CG
1966 11130 : max_iso_not0 = harmonics%max_iso_not0
1967 11130 : lmax_expansion = indso(1, max_iso_not0)
1968 11130 : max_s_harm = harmonics%max_s_harm
1969 :
1970 77910 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl))
1971 66780 : ALLOCATE (gVg_h(na, 0:2*maxl), gVg_s(na, 0:2*maxl))
1972 : ALLOCATE (matso_h(nsoset(maxl), nsoset(maxl)), &
1973 66780 : matso_s(nsoset(maxl), nsoset(maxl)))
1974 44520 : ALLOCATE (vx(na, nr))
1975 66780 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
1976 :
1977 713930 : g1 = 0.0_dp
1978 713930 : g2 = 0.0_dp
1979 11130 : m1 = 0
1980 37301 : DO iset1 = 1, nset
1981 26171 : n1 = nsoset(lmax(iset1))
1982 26171 : m2 = 0
1983 104104 : DO iset2 = 1, nset
1984 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
1985 77933 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
1986 77933 : CPASSERT(max_iso_not0_local .LE. max_iso_not0)
1987 :
1988 77933 : n2 = nsoset(lmax(iset2))
1989 285340 : DO ipgf1 = 1, npgf(iset1)
1990 207407 : ngau1 = n1*(ipgf1 - 1) + m1
1991 207407 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
1992 207407 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
1993 :
1994 12701757 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
1995 990514 : DO ipgf2 = 1, npgf(iset2)
1996 705174 : ngau2 = n2*(ipgf2 - 1) + m2
1997 :
1998 42454574 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
1999 705174 : lmin12 = lmin(iset1) + lmin(iset2)
2000 705174 : lmax12 = lmax(iset1) + lmax(iset2)
2001 :
2002 : ! reduce expansion local densities
2003 912581 : IF (lmin12 .LE. lmax_expansion) THEN
2004 :
2005 188648786 : gg = 0.0_dp
2006 704229 : IF (lmin12 == 0) THEN
2007 25158111 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2008 : ELSE
2009 17248268 : 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 704229 : IF (lmax12 .GT. lmax_expansion) lmax12 = lmax_expansion
2014 :
2015 989825 : DO l = lmin12 + 1, lmax12
2016 19770425 : gg(1:nr, l) = grid_atom%rad(1:nr)*gg(:, l - 1)
2017 : END DO
2018 :
2019 1563382 : DO ispin = 1, nspins
2020 859153 : ld = lmax12 + 1
2021 55319803 : DO ir = 1, nr
2022 2778352303 : 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 859153 : gg(1:nr, 0:lmax12), nr, 0.0_dp, gVg_h(1:na, 0:lmax12), na)
2026 55319803 : DO ir = 1, nr
2027 2778352303 : 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 859153 : gg(1:nr, 0:lmax12), nr, 0.0_dp, gVg_s(1:na, 0:lmax12), na)
2031 :
2032 80831721 : matso_h = 0.0_dp
2033 80831721 : matso_s = 0.0_dp
2034 6523692 : DO iso = 1, max_iso_not0_local
2035 16789867 : DO icg = 1, cg_n_list(iso)
2036 10266175 : iso1 = cg_list(1, icg, iso)
2037 10266175 : iso2 = cg_list(2, icg, iso)
2038 10266175 : l = indso(1, iso1) + indso(1, iso2)
2039 :
2040 10266175 : CPASSERT(l <= lmax_expansion)
2041 529239464 : 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 513308750 : 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 523574925 : harmonics%slm(ia, iso)
2050 : END DO
2051 : END DO
2052 : END DO
2053 :
2054 : ! Write in the global matrix
2055 3717179 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2056 2153797 : iso1 = nsoset(lmin(iset1) - 1) + 1
2057 2153797 : iso2 = ngau2 + ic
2058 : CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2059 2153797 : int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2060 : CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2061 3012950 : 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 182037 : m2 = m2 + maxso
2071 : END DO ! iset2
2072 37301 : m1 = m1 + maxso
2073 : END DO ! iset1
2074 :
2075 11130 : DEALLOCATE (g1, g2, gg, matso_h, matso_s, gVg_s, gVg_h, vx)
2076 :
2077 11130 : DEALLOCATE (cg_list, cg_n_list)
2078 :
2079 11130 : CALL timestop(handle)
2080 :
2081 11130 : 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 16648 : 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 16648 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2116 : REAL(dp) :: urad
2117 16648 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: g1, g2
2118 16648 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: dgg, gg, gVXCg_h, gVXCg_s, matso_h, &
2119 16648 : matso_s
2120 16648 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: gVXGg_h, gVXGg_s
2121 16648 : REAL(dp), DIMENSION(:, :), POINTER :: zet
2122 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
2123 : REAL(dp), DIMENSION(:, :, :, :), POINTER :: my_CG_dxyz
2124 :
2125 16648 : CALL timeset(routineN, handle)
2126 :
2127 16648 : 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 16648 : nset=nset, zet=zet)
2132 :
2133 16648 : nr = grid_atom%nr
2134 16648 : na = grid_atom%ng_sphere
2135 16648 : my_CG => harmonics%my_CG
2136 16648 : my_CG_dxyz => harmonics%my_CG_dxyz
2137 16648 : max_iso_not0 = harmonics%max_iso_not0
2138 16648 : lmax_expansion = indso(1, max_iso_not0)
2139 16648 : max_s_harm = harmonics%max_s_harm
2140 :
2141 149832 : ALLOCATE (g1(nr), g2(nr), gg(nr, 0:2*maxl), dgg(nr, 0:2*maxl))
2142 99888 : ALLOCATE (gVXCg_h(na, 0:2*maxl), gVXCg_s(na, 0:2*maxl))
2143 99888 : 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 149832 : 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 99888 : matso_s(nsoset(maxl), nsoset(maxl)))
2149 :
2150 35696 : DO ispin = 1, nspins
2151 :
2152 996528 : g1 = 0.0_dp
2153 996528 : g2 = 0.0_dp
2154 19048 : m1 = 0
2155 82715 : DO iset1 = 1, nset
2156 47019 : n1 = nsoset(lmax(iset1))
2157 47019 : m2 = 0
2158 206790 : DO iset2 = 1, nset
2159 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
2160 159771 : max_s_harm, lmax_expansion, cg_list, cg_n_list, max_iso_not0_local)
2161 159771 : 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 159771 : max_s_harm, lmax_expansion, dcg_list, dcg_n_list, dmax_iso_not0_local)
2164 :
2165 159771 : n2 = nsoset(lmax(iset2))
2166 498766 : DO ipgf1 = 1, npgf(iset1)
2167 338995 : ngau1 = n1*(ipgf1 - 1) + m1
2168 338995 : size1 = nsoset(lmax(iset1)) - nsoset(lmin(iset1) - 1)
2169 338995 : nngau1 = nsoset(lmin(iset1) - 1) + ngau1
2170 :
2171 18009825 : g1(1:nr) = EXP(-zet(ipgf1, iset1)*grid_atom%rad2(1:nr))
2172 1321824 : DO ipgf2 = 1, npgf(iset2)
2173 823058 : ngau2 = n2*(ipgf2 - 1) + m2
2174 :
2175 44002638 : g2(1:nr) = EXP(-zet(ipgf2, iset2)*grid_atom%rad2(1:nr))
2176 823058 : lmin12 = lmin(iset1) + lmin(iset2)
2177 823058 : lmax12 = lmax(iset1) + lmax(iset2)
2178 :
2179 : !test reduce expansion local densities
2180 823058 : IF (lmin12 .LE. lmax_expansion) THEN
2181 :
2182 185462284 : gg = 0.0_dp
2183 185462284 : dgg = 0.0_dp
2184 :
2185 822458 : IF (lmin12 == 0) THEN
2186 30060597 : gg(1:nr, lmin12) = g1(1:nr)*g2(1:nr)
2187 : ELSE
2188 13911441 : 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 822458 : IF (lmax12 .GT. lmax_expansion) lmax12 = lmax_expansion
2193 :
2194 1301010 : DO l = lmin12 + 1, lmax12
2195 25864952 : 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 26687410 : 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 43972038 : gg(1:nr, lmax12)
2202 :
2203 175698560 : gVXCg_h = 0.0_dp
2204 175698560 : gVXCg_s = 0.0_dp
2205 690031268 : gVXGg_h = 0.0_dp
2206 690031268 : gVXGg_s = 0.0_dp
2207 :
2208 : ! Cross Term
2209 2123468 : DO l = lmin12, lmax12
2210 67125200 : DO ia = 1, na
2211 3496644142 : 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 3430341400 : 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 3430341400 : vxg_s(3, ia, ir, ispin)*harmonics%a(3, ia))
2225 :
2226 3430341400 : 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 3430341400 : gg(ir, l)*urad
2231 :
2232 : gVXGg_h(2, ia, l) = gVXGg_h(2, ia, l) + &
2233 : vxg_h(2, ia, ir, ispin)* &
2234 3430341400 : gg(ir, l)*urad
2235 :
2236 : gVXGg_h(3, ia, l) = gVXGg_h(3, ia, l) + &
2237 : vxg_h(3, ia, ir, ispin)* &
2238 3430341400 : gg(ir, l)*urad
2239 :
2240 : gVXGg_s(1, ia, l) = gVXGg_s(1, ia, l) + &
2241 : vxg_s(1, ia, ir, ispin)* &
2242 3430341400 : gg(ir, l)*urad
2243 :
2244 : gVXGg_s(2, ia, l) = gVXGg_s(2, ia, l) + &
2245 : vxg_s(2, ia, ir, ispin)* &
2246 3430341400 : gg(ir, l)*urad
2247 :
2248 : gVXGg_s(3, ia, l) = gVXGg_s(3, ia, l) + &
2249 : vxg_s(3, ia, ir, ispin)* &
2250 3495343132 : gg(ir, l)*urad
2251 :
2252 : END DO ! ir
2253 : END DO ! ia
2254 : END DO ! l
2255 :
2256 65854302 : matso_h = 0.0_dp
2257 65854302 : matso_s = 0.0_dp
2258 5476224 : DO iso = 1, max_iso_not0_local
2259 14629934 : DO icg = 1, cg_n_list(iso)
2260 9153710 : iso1 = cg_list(1, icg, iso)
2261 9153710 : iso2 = cg_list(2, icg, iso)
2262 :
2263 9153710 : l = indso(1, iso1) + indso(1, iso2)
2264 :
2265 : !test reduce expansion local densities
2266 9153710 : CPASSERT(l <= lmax_expansion)
2267 471244636 : DO ia = 1, na
2268 : matso_h(iso1, iso2) = matso_h(iso1, iso2) + &
2269 : gVXCg_h(ia, l)* &
2270 : harmonics%slm(ia, iso)* &
2271 457437160 : 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 466590870 : 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 2913706 : DO iso = 1, dmax_iso_not0_local
2285 16847874 : DO icg = 1, dcg_n_list(iso)
2286 13934168 : iso1 = dcg_list(1, icg, iso)
2287 13934168 : iso2 = dcg_list(2, icg, iso)
2288 :
2289 13934168 : l = indso(1, iso1) + indso(1, iso2)
2290 : !test reduce expansion local densities
2291 13934168 : CPASSERT(l <= lmax_expansion)
2292 711994484 : 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 695969068 : 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 709903236 : 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 3167403 : DO ic = nsoset(lmin(iset2) - 1) + 1, nsoset(lmax(iset2))
2316 2005350 : iso1 = nsoset(lmin(iset1) - 1) + 1
2317 2005350 : iso2 = ngau2 + ic
2318 : CALL daxpy(size1, 1.0_dp, matso_h(iso1, ic), 1, &
2319 2005350 : int_hh(ispin)%r_coef(nngau1 + 1, iso2), 1)
2320 : CALL daxpy(size1, 1.0_dp, matso_s(iso1, ic), 1, &
2321 2828408 : int_ss(ispin)%r_coef(nngau1 + 1, iso2), 1)
2322 : END DO
2323 :
2324 : END DO ! ipfg2
2325 : END DO ! ipfg1
2326 526332 : m2 = m2 + maxso
2327 : END DO ! iset2
2328 66067 : m1 = m1 + maxso
2329 : END DO ! iset1
2330 : END DO ! ispin
2331 :
2332 16648 : DEALLOCATE (g1, g2, gg, dgg, matso_h, matso_s, gVXCg_h, gVXCg_s, gVXGg_h, gVXGg_s)
2333 16648 : DEALLOCATE (cg_list, cg_n_list, dcg_list, dcg_n_list)
2334 :
2335 16648 : CALL timestop(handle)
2336 :
2337 16648 : 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 308 : 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 308 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf
2368 308 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: fga, fgr, r1, r2, work
2369 308 : REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: a1, a2, intso_h, intso_s
2370 308 : REAL(dp), DIMENSION(:, :), POINTER :: slm, zet
2371 308 : REAL(dp), DIMENSION(:, :, :), POINTER :: dslm_dxyz
2372 :
2373 308 : CALL timeset(routineN, handle)
2374 :
2375 308 : 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 308 : maxso=maxso, npgf=npgf, nset=nset, zet=zet)
2379 :
2380 308 : na = grid_atom%ng_sphere
2381 308 : nr = grid_atom%nr
2382 :
2383 308 : slm => harmonics%slm
2384 308 : 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 2156 : ALLOCATE (a1(na, nset*maxso, 3), a2(na, nset*maxso, 3))
2389 1848 : ALLOCATE (r1(nr, nset*maxso), r2(nr, nset*maxso))
2390 3859916 : a1 = 0.0_dp; a2 = 0.0_dp
2391 1300788 : r1 = 0.0_dp; r2 = 0.0_dp
2392 :
2393 1008 : DO iset = 1, nset
2394 3200 : DO ipgf = 1, npgf(iset)
2395 2192 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
2396 8211 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
2397 5319 : 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 289469 : 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 289469 : *EXP(-zet(ipgf, iset)*grid_atom%rad2(1:nr))
2408 :
2409 23468 : DO dir = 1, 3
2410 : ! the purely angular part of d/dx(r^l Y_lm) * exp(-al*r^2)
2411 853119 : 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 858438 : 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 1540 : ALLOCATE (intso_h(nset*maxso, nset*maxso, nspins))
2423 1232 : ALLOCATE (intso_s(nset*maxso, nset*maxso, nspins))
2424 5111456 : intso_h = 0.0_dp; intso_s = 0.0_dp
2425 :
2426 924 : ALLOCATE (fga(na, 1))
2427 924 : ALLOCATE (fgr(nr, 1))
2428 616 : ALLOCATE (work(na, 1))
2429 50604 : fga = 0.0_dp; fgr = 0.0_dp; work = 0.0_dp
2430 :
2431 1008 : DO iset = 1, nset
2432 2984 : DO jset = 1, nset
2433 9254 : DO ipgf = 1, npgf(iset)
2434 6578 : starti = (iset - 1)*maxso + (ipgf - 1)*nsoset(lmax(iset))
2435 33710 : DO jpgf = 1, npgf(jset)
2436 25156 : startj = (jset - 1)*maxso + (jpgf - 1)*nsoset(lmax(jset))
2437 :
2438 99450 : DO iso = nsoset(lmin(iset) - 1) + 1, nsoset(lmax(iset))
2439 296525 : 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 10704803 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r1(1:nr, startj + jso)
2445 814612 : DO dir = 1, 3
2446 31846869 : fga(1:na, 1) = a1(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
2447 :
2448 1425571 : DO ispin = 1, nspins
2449 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2450 610959 : 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 610959 : 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 610959 : 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 1221918 : 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 10704803 : fgr(1:nr, 1) = r1(1:nr, starti + iso)*r2(1:nr, startj + jso)
2463 814612 : DO dir = 1, 3
2464 31846869 : fga(1:na, 1) = a1(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
2465 :
2466 1425571 : DO ispin = 1, nspins
2467 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2468 610959 : 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 610959 : 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 610959 : 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 1221918 : 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 10704803 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r1(1:nr, startj + jso)
2481 814612 : DO dir = 1, 3
2482 31846869 : fga(1:na, 1) = a2(1:na, starti + iso, dir)*a1(1:na, startj + jso, dir)
2483 :
2484 1425571 : DO ispin = 1, nspins
2485 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2486 610959 : 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 610959 : 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 610959 : 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 1221918 : 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 10704803 : fgr(1:nr, 1) = r2(1:nr, starti + iso)*r2(1:nr, startj + jso)
2499 882328 : DO dir = 1, 3
2500 31846869 : fga(1:na, 1) = a2(1:na, starti + iso, dir)*a2(1:na, startj + jso, dir)
2501 :
2502 1425571 : DO ispin = 1, nspins
2503 : CALL dgemm('N', 'N', na, 1, nr, 0.5_dp, vtau_h(:, :, ispin), na, fgr, &
2504 610959 : 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 610959 : 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 610959 : 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 1221918 : 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 616 : DO ispin = 1, nspins
2525 2555574 : int_hh(ispin)%r_coef(:, :) = int_hh(ispin)%r_coef(:, :) + intso_h(:, :, ispin)
2526 2555882 : int_ss(ispin)%r_coef(:, :) = int_ss(ispin)%r_coef(:, :) + intso_s(:, :, ispin)
2527 : END DO
2528 :
2529 308 : CALL timestop(handle)
2530 :
2531 616 : END SUBROUTINE dgaVtaudgb
2532 :
2533 : END MODULE qs_vxc_atom
|
