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