(git:40b9a06)
Loading...
Searching...
No Matches
qs_vxc_atom.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief routines that build the integrals of the Vxc potential calculated
10!> for the atomic density in the basis set of spherical primitives
11! **************************************************************************************************
12MODULE qs_vxc_atom
13 USE atomic_kind_types, ONLY: atomic_kind_type,&
14 get_atomic_kind
15 USE basis_set_types, ONLY: gto_basis_set_type
16 USE cell_types, ONLY: cell_type
17 USE cp_array_utils, ONLY: cp_3d_r_cp_type
18 USE cp_control_types, ONLY: dft_control_type
19 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
20 USE external_potential_types, ONLY: gth_potential_type,&
21 sgp_potential_type
22 USE input_constants, ONLY: cdft_alpha_constraint,&
23 cdft_beta_constraint,&
24 cdft_charge_constraint,&
25 cdft_magnetization_constraint,&
26 xc_none
27 USE input_section_types, ONLY: section_get_rval,&
28 section_vals_get_subs_vals,&
29 section_vals_type,&
30 section_vals_val_get
31 USE kinds, ONLY: dp,&
32 int_8
33 USE lebedev, ONLY: deallocate_lebedev_grids,&
34 get_number_of_lebedev_grid,&
35 init_lebedev_grids,&
36 lebedev_grid
37 USE memory_utilities, ONLY: reallocate
38 USE message_passing, ONLY: mp_para_env_type
39 USE orbital_pointers, ONLY: indco,&
40 indso,&
41 nco,&
42 ncoset,&
43 nsoset
44 USE orbital_transformation_matrices, ONLY: orbtramat
45 USE particle_types, ONLY: particle_type
46 USE pw_env_types, ONLY: pw_env_get,&
47 pw_env_type
48 USE pw_grid_types, ONLY: pw_grid_type
49 USE pw_methods, ONLY: pw_axpy
50 USE pw_pool_types, ONLY: pw_pool_type
51 USE pw_types, ONLY: pw_c1d_gs_type,&
52 pw_r3d_rs_type
53 USE qs_cdft_grid, ONLY: cdft_point_context_create,&
54 cdft_point_context_release,&
55 cdft_point_context_type,&
56 cdft_point_weights
57 USE qs_cdft_types, ONLY: cdft_control_type
58 USE qs_environment_types, ONLY: get_qs_env,&
59 qs_environment_type
60 USE qs_force_types, ONLY: qs_force_type
61 USE qs_grid_atom, ONLY: allocate_grid_atom,&
62 create_grid_atom,&
63 grid_atom_type
64 USE qs_harmonics_atom, ONLY: harmonics_atom_type
65 USE qs_kind_types, ONLY: get_qs_kind,&
66 has_nlcc,&
67 qs_kind_type
68 USE qs_linres_types, ONLY: nablavks_atom_type
69 USE qs_native_grid_cache, ONLY: fetch_native_grid_stencil,&
70 native_grid_interp_npts,&
71 native_grid_interp_offset_max,&
72 native_grid_interp_offset_min,&
73 native_grid_interpolation_stencil_type,&
74 prepare_native_grid_cache,&
75 store_native_grid_stencil
76 USE qs_rho_atom_methods, ONLY: replicate_rho_atom_radial
77 USE qs_rho_atom_types, ONLY: get_rho_atom,&
78 rho_atom_coeff,&
79 rho_atom_type
80 USE qs_rho_types, ONLY: qs_rho_get,&
81 qs_rho_type
82 USE qs_vxc_atom_utils, ONLY: &
83 calc_rho_angular, calc_rho_nlcc, calc_tau_atom, calc_weight_function, &
84 create_tau_basis_cache, dgavtaudgb, evaluate_nlcc_primitive_fields, gavxcgb_gc, &
85 gavxcgb_nogc, gapw_atom_grid_support_radius, release_tau_basis_cache, tau_basis_cache_type
86 USE skala_gpw_features, ONLY: build_periodic_atom_image_layout,&
87 periodic_atom_image_partition,&
88 periodic_atom_image_partition_from_layout,&
89 skala_gpw_smooth_partition_derivatives,&
90 smooth_atom_partition,&
91 smooth_partition_atomic_weight_scale,&
92 smooth_partition_atomic_weight_scale_derivative
93 USE skala_gpw_functional, ONLY: &
94 build_vxc_from_feature_grads, get_gauxc_section, native_skala_gapw_density_partition, &
95 skala_gapw_atom_composite_energy, skala_gapw_atom_vxc_of_r, &
96 skala_gapw_density_partition_hard_minus_soft, skala_gapw_density_partition_hard_only, &
97 skala_gapw_density_partition_none, skala_gapw_density_partition_soft_only, &
98 xc_section_uses_gauxc_model, xc_section_uses_native_skala_evaluator
99 USE spherical_harmonics, ONLY: y_lm
100 USE util, ONLY: get_limit
101 USE virial_types, ONLY: virial_type
102 USE xc_atom, ONLY: fill_rho_set,&
103 vxc_of_r_epr,&
104 vxc_of_r_new,&
105 xc_rho_set_atom_update
106 USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
107 xc_dset_create,&
108 xc_dset_release,&
109 xc_dset_zero_all
110 USE xc_derivatives, ONLY: xc_functionals_get_needs
111 USE xc_input_constants, ONLY: skala_gapw_cp2k_default,&
112 skala_gapw_direct_valence,&
113 skala_gapw_paw_one_center,&
114 skala_gapw_paw_one_center_split
115 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
116 USE xc_rho_set_types, ONLY: xc_rho_set_create,&
117 xc_rho_set_get,&
118 xc_rho_set_release,&
119 xc_rho_set_type,&
120 xc_rho_set_update
121#include "./base/base_uses.f90"
122
123 IMPLICIT NONE
124
125 PRIVATE
126
127 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_vxc_atom'
128
129 ! A wrapped stencil can touch both end tiles and one adjacent interior tile when the final
130 ! tile is shorter than the stencil. Three tiles per direction are therefore sufficient.
131 INTEGER, PARAMETER, PRIVATE :: native_grid_adjoint_tile_edge = 64, &
132 native_grid_adjoint_max_tiles_per_direction = 3, &
133 native_grid_adjoint_max_bins_per_row = &
134 native_grid_adjoint_max_tiles_per_direction**3
135
136 PUBLIC :: calculate_vxc_atom, &
137 calculate_vxc_atom_epr, &
138 gapw_cdft_one_center
139
140! **************************************************************************************************
141
142CONTAINS
143
144! **************************************************************************************************
145!> \brief Initialize atom-centered quadrature for native Skala layouts.
146!> \param kind_set quantum kinds
147!> \param dft_control DFT controls supplying the radial quadrature
148! **************************************************************************************************
149 SUBROUTINE ensure_native_skala_atom_grids(kind_set, dft_control)
150 TYPE(qs_kind_type), DIMENSION(:), POINTER :: kind_set
151 TYPE(dft_control_type), POINTER :: dft_control
152
153 INTEGER :: ikind, ll, na, nr, quadrature
154 TYPE(grid_atom_type), POINTER :: grid_atom
155
156 quadrature = dft_control%qs_control%gapw_control%quadrature
157 CALL init_lebedev_grids()
158 DO ikind = 1, SIZE(kind_set)
159 NULLIFY (grid_atom)
160 CALL get_qs_kind(kind_set(ikind), grid_atom=grid_atom, ngrid_ang=na, ngrid_rad=nr)
161 IF (ASSOCIATED(grid_atom)) THEN
162 IF (ASSOCIATED(grid_atom%weight) .AND. grid_atom%nr == nr) cycle
163 ELSE
164 CALL allocate_grid_atom(kind_set(ikind)%grid_atom)
165 grid_atom => kind_set(ikind)%grid_atom
166 END IF
167 ll = get_number_of_lebedev_grid(n=na)
168 na = lebedev_grid(ll)%n
169 grid_atom%ng_sphere = na
170 grid_atom%nr = nr
171 CALL create_grid_atom(grid_atom, nr, na, 0, ll, quadrature)
172 END DO
173 CALL deallocate_lebedev_grids()
174
175 END SUBROUTINE ensure_native_skala_atom_grids
176
177! **************************************************************************************************
178!> \brief Decide whether a kind contributes hard-minus-soft primitive fields.
179!> \param paw_atom whether CP2K constructed a one-center representation for the kind
180!> \param gapw_representation requested Skala pseudopotential GAPW representation
181!> \param has_pseudopotential whether the kind uses a GTH or semi-global pseudopotential
182!> \param zeff valence charge of the potential
183!> \param zatom atomic number
184!> \return true when hard-minus-soft fields contribute for this kind
185! **************************************************************************************************
186 PURE FUNCTION native_skala_uses_one_center_kind( &
187 paw_atom, gapw_representation, has_pseudopotential, zeff, zatom) RESULT(use_one_center)
188 LOGICAL, INTENT(IN) :: paw_atom
189 INTEGER, INTENT(IN) :: gapw_representation
190 LOGICAL, INTENT(IN) :: has_pseudopotential
191 REAL(dp), INTENT(IN) :: zeff
192 INTEGER, INTENT(IN) :: zatom
193 LOGICAL :: use_one_center
194
195 use_one_center = paw_atom
196 IF (.NOT. use_one_center) RETURN
197
198 SELECT CASE (gapw_representation)
199 CASE (skala_gapw_direct_valence)
200 use_one_center = .NOT. has_pseudopotential
201 CASE (skala_gapw_paw_one_center, skala_gapw_paw_one_center_split)
202 CONTINUE
203 CASE (skala_gapw_cp2k_default)
204 IF (has_pseudopotential .AND. &
205 abs(zeff - real(zatom, dp)) <= 1.0e-10_dp) use_one_center = .false.
206 END SELECT
207
208 END FUNCTION native_skala_uses_one_center_kind
209
210! **************************************************************************************************
211!> \brief ...
212!> \param qs_env ...
213!> \param energy_only ...
214!> \param exc1 the on-body ex energy contribution
215!> \param adiabatic_rescale_factor ...
216!> \param kind_set_external provides a non-default kind_set to use
217!> \param rho_atom_set_external provides a non-default atomic density set to use
218!> \param xc_section_external provides an external non-default XC
219!> \param calculate_forces ...
220!> \param composite_vxc_rho ...
221!> \param composite_vxc_tau ...
222!> \param composite_reference_active ...
223!> \param direct_valence_atom_grid evaluate the smooth valence fields on atom-centered grids
224!> \param atom_composite_grid evaluate GAPW primitive fields on atom-centered composite grids
225! **************************************************************************************************
226 SUBROUTINE calculate_vxc_atom(qs_env, energy_only, exc1, &
227 adiabatic_rescale_factor, kind_set_external, &
228 rho_atom_set_external, xc_section_external, calculate_forces, &
229 composite_vxc_rho, composite_vxc_tau, composite_reference_active, &
230 direct_valence_atom_grid, atom_composite_grid)
231
232 TYPE(qs_environment_type), POINTER :: qs_env
233 LOGICAL, INTENT(IN) :: energy_only
234 REAL(dp), INTENT(INOUT) :: exc1
235 REAL(dp), INTENT(IN), OPTIONAL :: adiabatic_rescale_factor
236 TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
237 POINTER :: kind_set_external
238 TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
239 POINTER :: rho_atom_set_external
240 TYPE(section_vals_type), OPTIONAL, POINTER :: xc_section_external
241 LOGICAL, INTENT(IN), OPTIONAL :: calculate_forces
242 TYPE(pw_r3d_rs_type), DIMENSION(:), OPTIONAL, &
243 POINTER :: composite_vxc_rho, composite_vxc_tau
244 LOGICAL, INTENT(OUT), OPTIONAL :: composite_reference_active
245 LOGICAL, INTENT(IN), OPTIONAL :: direct_valence_atom_grid, &
246 atom_composite_grid
247
248 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_vxc_atom'
249
250 INTEGER :: adjoint_bin, adjoint_entry, adjoint_nbins, adjoint_nchannels, &
251 adjoint_tile_count(3), adjoint_tile_lower(3), adjoint_tile_upper(3), &
252 atom_composite_components, bo(2), composite_descriptor_target_image, &
253 composite_image_periodicity(3), composite_local_atom, composite_local_natom, &
254 composite_nflat, composite_partition_target_image, composite_pw_nflat, composite_row, &
255 gapw_density_partition, gapw_representation, handle, ia, iat, iatom, icomponent, idir, &
256 ikind, image_i1, image_i2, image_i3, image_lower(3), image_shell(3), image_shift(3), &
257 image_upper(3), ir, ispin, iw, jdir, myfun, na
258 INTEGER :: natom, nr, nspins, num_pe, on, source_atom, target_atom, xc_deriv_method_id, &
259 xc_rho_smooth_id, zatom
260 INTEGER(KIND=int_8), ALLOCATABLE, DIMENSION(:) :: composite_atomic_grid_sizes, &
261 composite_local_grid_sizes
262 INTEGER, ALLOCATABLE, DIMENSION(:) :: adjoint_bin_offsets, adjoint_bin_rows, &
263 composite_atom_end, composite_atom_kind, composite_atom_kind_index, composite_atom_start, &
264 composite_grid_atom, composite_local_atoms
265 INTEGER, DIMENSION(2, 3) :: bounds
266 INTEGER, DIMENSION(:), POINTER :: atom_list
267 LOGICAL :: accint, atom_composite_active, atom_composite_diagnostic, &
268 atom_composite_reference, direct_valence_atom_composite, donlcc, evaluate_hard, &
269 evaluate_soft, gradient_f, image_partition_atom_composite, lsd, my_calculate_forces, &
270 native_grid_diagnostics, nlcc, one_center_kind, paw_atom, paw_pseudopotentials, &
271 requested_atom_composite_grid, rho_g_valid, skala_atom_grid, source_matrix_local, tau_f, &
272 tau_r_valid, use_atom_composite_density, use_atom_composite_gradient, &
273 use_atom_composite_tau, use_virial
274 LOGICAL, ALLOCATABLE, DIMENSION(:) :: composite_partition_included
275 REAL(dp) :: agr, alpha, atom_composite_exc, atom_composite_nelec, &
276 composite_cross_cutoff_max, composite_cross_density_max, composite_cross_grad_max, &
277 composite_cross_kin_max, composite_density_max, composite_density_min, &
278 composite_grad_max, composite_kin_max, composite_kin_min, composite_tau_integral, &
279 cross_cutoff, density_cut, descriptor_window_adjoint, descriptor_window_weight, exc_h, &
280 exc_s, feature_vxc_analytic, feature_vxc_fd, feature_vxc_minus, feature_vxc_plus, &
281 feature_vxc_step, gradient_cut, local_partition_weight, my_adiabatic_rescale_factor, &
282 nlcc_density, nlcc_spin_factor
283 REAL(dp) :: one_center_density_field_contraction, one_center_density_matrix_contraction, &
284 one_center_field_contraction, one_center_gradient_field_contraction, &
285 one_center_gradient_matrix_contraction, one_center_matrix_contraction, &
286 one_center_rho_grad_field_contraction, one_center_rho_grad_matrix_contraction, &
287 one_center_tau_field_contraction, one_center_tau_matrix_contraction, &
288 one_center_tensor_contraction, partition_adjoint, partition_scale, partition_weight, &
289 smooth_grid_contraction, smooth_input_contraction, target_partition_adjoint, tau_cut, zeff
290 REAL(dp), ALLOCATABLE, DIMENSION(:) :: composite_atomic_grid_weight_grad, &
291 composite_atomic_grid_weights, composite_base_grid_weights, composite_distances, &
292 composite_grid_weight_grad, composite_grid_weights, composite_partition_weights, fw
293 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: composite_atom_coord_grad, composite_atom_coords, &
294 composite_cross_density, composite_cross_force, composite_cross_force_local, &
295 composite_cross_kin, composite_density, composite_density_grad, &
296 composite_descriptor_image_coords, composite_explicit_force, composite_grid_coord_force, &
297 composite_grid_coord_grad, composite_grid_coords, composite_kin, composite_kin_grad, &
298 composite_local_atom_coords, composite_model_atom_force, composite_moving_smooth_force, &
299 composite_nlcc_center_force, composite_nlcc_center_force_local, &
300 composite_nlcc_target_force
301 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: composite_partition_atom_coords, &
302 composite_partition_force, composite_partition_force_local, &
303 composite_partition_image_coords, composite_smooth_density_cache, &
304 composite_smooth_kin_cache, local_partition_datom
305 REAL(dp), ALLOCATABLE, DIMENSION(:, :), TARGET :: smooth_density_adjoint_storage, &
306 smooth_kin_adjoint_storage
307 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :) :: composite_cross_grad, composite_grad, &
308 composite_grad_grad, composite_int_h, composite_int_s, composite_partition_datom, &
309 composite_partition_dstrain, composite_smooth_gradient_cache, vtau_h_local, vtau_s_local, &
310 vxc_h_local, vxc_s_local
311 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :), TARGET :: smooth_grad_adjoint_storage
312 REAL(dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: vxg_h_local, vxg_s_local
313 REAL(dp), DIMENSION(1, 1, 1) :: tau_d
314 REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
315 REAL(dp), DIMENSION(2) :: composite_smooth_density_adjoint_value, &
316 composite_smooth_density_value, composite_smooth_kin_adjoint_value, &
317 composite_smooth_kin_value, cross_density, cross_density_adjoint, cross_kin, &
318 cross_kin_adjoint
319 REAL(dp), DIMENSION(3) :: composite_point, cross_displacement, cross_spatial_derivative, &
320 fractional, image_translation, nlcc_gradient, nlcc_spatial_derivative, &
321 skala_atom_force_h, skala_atom_force_s, spatial_derivative
322 REAL(dp), DIMENSION(3, 1) :: local_descriptor_datom
323 REAL(dp), DIMENSION(3, 2) :: composite_smooth_gradient_adjoint_value, &
324 composite_smooth_gradient_value, cross_density_spatial, cross_grad, cross_grad_adjoint, &
325 cross_kin_spatial
326 REAL(dp), DIMENSION(3, 3) :: composite_cross_image_virial, &
327 composite_cross_image_virial_local, composite_explicit_virial, composite_feature_virial, &
328 composite_interpolation_virial, composite_partition_strain_virial, &
329 local_descriptor_dstrain, local_partition_dstrain, nlcc_hessian, skala_atom_virial, &
330 skala_atom_virial_h, skala_atom_virial_s
331 REAL(dp), DIMENSION(3, 3, 2) :: cross_grad_spatial
332 REAL(dp), DIMENSION(4) :: feature_component_analytic, &
333 feature_component_fd
334 REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, smooth_density_adjoint, &
335 smooth_kin_adjoint, weight_h, weight_s
336 REAL(dp), DIMENSION(:, :, :), POINTER :: composite_smooth_rho, composite_smooth_rhoa, &
337 composite_smooth_rhob, composite_smooth_tau, composite_smooth_tau_a, &
338 composite_smooth_tau_b, rho_h, rho_s, smooth_grad_adjoint, smooth_rho, smooth_rhoa, &
339 smooth_rhob, smooth_tau, smooth_tau_a, smooth_tau_b, tau_h, tau_s, vtau_h, vtau_s, vxc_h, &
340 vxc_s
341 REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
342 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
343 TYPE(cell_type), POINTER :: cell
344 TYPE(cp_3d_r_cp_type), DIMENSION(3) :: composite_smooth_drho, composite_smooth_drhoa, &
345 composite_smooth_drhob, smooth_drho, smooth_drhoa, smooth_drhob
346 TYPE(dft_control_type), POINTER :: dft_control
347 TYPE(grid_atom_type), POINTER :: grid_atom
348 TYPE(gth_potential_type), POINTER :: gth_potential
349 TYPE(gto_basis_set_type), POINTER :: basis_1c
350 TYPE(harmonics_atom_type), POINTER :: harmonics
351 TYPE(mp_para_env_type), POINTER :: para_env
352 TYPE(native_grid_interpolation_stencil_type) :: interpolation_stencil
353 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
354 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: smooth_rho_g
355 TYPE(pw_env_type), POINTER :: pw_env
356 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
357 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: smooth_rho_r, smooth_tau_r, &
358 smooth_vxc_rho, smooth_vxc_tau
359 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
360 TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
361 TYPE(qs_rho_type), POINTER :: rho_struct
362 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s, dr_h, dr_s, int_hh, &
363 int_ss, r_h, r_s
364 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
365 TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom_set
366 TYPE(rho_atom_type), POINTER :: rho_atom
367 TYPE(section_vals_type), POINTER :: gauxc_section, input, my_xc_section, &
368 xc_fun_section
369 TYPE(sgp_potential_type), POINTER :: sgp_potential
370 TYPE(tau_basis_cache_type) :: tau_basis_cache
371 TYPE(virial_type), POINTER :: virial
372 TYPE(xc_derivative_set_type) :: deriv_set
373 TYPE(xc_rho_cflags_type) :: needs
374 TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s, smooth_rho_set
375
376! -------------------------------------------------------------------------
377
378 CALL timeset(routinen, handle)
379
380 NULLIFY (atom_list)
381 NULLIFY (auxbas_pw_pool)
382 NULLIFY (my_kind_set)
383 NULLIFY (atomic_kind_set)
384 NULLIFY (cell)
385 NULLIFY (grid_atom)
386 NULLIFY (gth_potential)
387 NULLIFY (force)
388 NULLIFY (harmonics)
389 NULLIFY (input)
390 NULLIFY (para_env)
391 NULLIFY (particle_set)
392 NULLIFY (pw_env)
393 NULLIFY (rho_atom)
394 NULLIFY (rho_struct)
395 NULLIFY (my_rho_atom_set)
396 NULLIFY (rho_nlcc)
397 NULLIFY (smooth_rho, smooth_rhoa, smooth_rhob, smooth_tau, smooth_tau_a, smooth_tau_b)
398 NULLIFY (composite_smooth_rho, composite_smooth_rhoa, composite_smooth_rhob, &
399 composite_smooth_tau, composite_smooth_tau_a, composite_smooth_tau_b)
400 NULLIFY (smooth_rho_g, smooth_rho_r, smooth_tau_r)
401 NULLIFY (smooth_vxc_rho, smooth_vxc_tau)
402 DO idir = 1, 3
403 NULLIFY (smooth_drho(idir)%array, smooth_drhoa(idir)%array, smooth_drhob(idir)%array)
404 NULLIFY (composite_smooth_drho(idir)%array, &
405 composite_smooth_drhoa(idir)%array, &
406 composite_smooth_drhob(idir)%array)
407 END DO
408 NULLIFY (sgp_potential)
409 NULLIFY (virial)
410 my_calculate_forces = .false.
411 IF (PRESENT(calculate_forces)) my_calculate_forces = calculate_forces
412 IF (PRESENT(composite_reference_active)) composite_reference_active = .false.
413 direct_valence_atom_composite = .false.
414 IF (PRESENT(direct_valence_atom_grid)) THEN
415 direct_valence_atom_composite = direct_valence_atom_grid
416 END IF
417 requested_atom_composite_grid = .false.
418 IF (PRESENT(atom_composite_grid)) requested_atom_composite_grid = atom_composite_grid
419
420 IF (PRESENT(adiabatic_rescale_factor)) THEN
421 my_adiabatic_rescale_factor = adiabatic_rescale_factor
422 ELSE
423 my_adiabatic_rescale_factor = 1.0_dp
424 END IF
425
426 CALL get_qs_env(qs_env=qs_env, &
427 dft_control=dft_control, &
428 cell=cell, &
429 para_env=para_env, &
430 atomic_kind_set=atomic_kind_set, &
431 qs_kind_set=my_kind_set, &
432 input=input, &
433 particle_set=particle_set, &
434 pw_env=pw_env, &
435 virial=virial, &
436 rho_atom_set=my_rho_atom_set, &
437 force=force)
438
439 IF (dft_control%qs_control%gapw_xc) THEN
440 CALL get_qs_env(qs_env=qs_env, rho_xc=rho_struct)
441 ELSE
442 CALL get_qs_env(qs_env=qs_env, rho=rho_struct)
443 END IF
444
445 IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
446 IF (PRESENT(rho_atom_set_external)) my_rho_atom_set => rho_atom_set_external
447
448 nlcc = has_nlcc(my_kind_set)
449 accint = dft_control%qs_control%gapw_control%accurate_xcint
450
451 my_xc_section => section_vals_get_subs_vals(input, "DFT%XC")
452
453 IF (PRESENT(xc_section_external)) my_xc_section => xc_section_external
454
455 xc_fun_section => section_vals_get_subs_vals(my_xc_section, "XC_FUNCTIONAL")
456 CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", &
457 i_val=myfun)
458 skala_atom_grid = xc_section_uses_gauxc_model(my_xc_section)
459 gapw_representation = skala_gapw_cp2k_default
460 atom_composite_diagnostic = .false.
461 atom_composite_reference = .false.
462 paw_pseudopotentials = .false.
463 native_grid_diagnostics = .false.
464 atom_composite_components = 1
465 feature_vxc_step = 3.0e-3_dp
466 IF (skala_atom_grid) THEN
467 gauxc_section => get_gauxc_section(my_xc_section)
468 cpassert(ASSOCIATED(gauxc_section))
469 CALL section_vals_val_get(gauxc_section, "PSEUDOPOTENTIAL_GAPW_REPRESENTATION", &
470 i_val=gapw_representation)
471 CALL section_vals_val_get(gauxc_section, "NATIVE_GRID_DIAGNOSTICS", &
472 l_val=native_grid_diagnostics)
473 END IF
474 IF (skala_atom_grid) THEN
475 DO ikind = 1, SIZE(my_kind_set)
476 NULLIFY (gth_potential, sgp_potential)
477 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
478 gth_potential=gth_potential, sgp_potential=sgp_potential)
479 paw_pseudopotentials = paw_pseudopotentials .OR. &
480 (paw_atom .AND. (ASSOCIATED(gth_potential) .OR. &
481 ASSOCIATED(sgp_potential)))
482 END DO
483 END IF
484 IF (skala_atom_grid .AND. xc_section_uses_native_skala_evaluator(my_xc_section)) THEN
485 CALL section_vals_val_get(gauxc_section, &
486 "NATIVE_GRID_GAPW_ATOM_COMPOSITE_DIAGNOSTIC", &
487 l_val=atom_composite_diagnostic)
488 CALL section_vals_val_get(gauxc_section, &
489 "NATIVE_GRID_GAPW_ATOM_COMPOSITE_REFERENCE", &
490 l_val=atom_composite_reference)
491 CALL section_vals_val_get(gauxc_section, &
492 "NATIVE_GRID_GAPW_ATOM_COMPOSITE_COMPONENTS", &
493 i_val=atom_composite_components)
494 CALL section_vals_val_get(gauxc_section, &
495 "NATIVE_GRID_GAPW_ATOM_COMPOSITE_FD_STEP", &
496 r_val=feature_vxc_step)
497 END IF
498 atom_composite_reference = atom_composite_reference .OR. &
499 (gapw_representation == skala_gapw_paw_one_center .AND. &
500 paw_pseudopotentials)
501 atom_composite_reference = atom_composite_reference .OR. requested_atom_composite_grid
502 atom_composite_reference = atom_composite_reference .OR. direct_valence_atom_composite
503 atom_composite_active = atom_composite_diagnostic .OR. atom_composite_reference
504 use_atom_composite_density = atom_composite_components <= 2
505 use_atom_composite_gradient = atom_composite_components <= 2
506 use_atom_composite_tau = atom_composite_components == 1 .OR. &
507 atom_composite_components == 3
508 IF (atom_composite_active) THEN
509 CALL ensure_native_skala_atom_grids(my_kind_set, dft_control)
510 END IF
511 IF (direct_valence_atom_composite) THEN
512 use_atom_composite_density = .false.
513 use_atom_composite_gradient = .false.
514 use_atom_composite_tau = .false.
515 END IF
516 ! CP2K's auxiliary PW fields are represented on one index-periodic cell even
517 ! when the physical Poisson problem is isolated or partially periodic.
518 image_partition_atom_composite = atom_composite_active
519 composite_image_periodicity = 1
520 IF (PRESENT(composite_reference_active)) composite_reference_active = atom_composite_reference
521 gapw_density_partition = skala_gapw_density_partition_hard_minus_soft
522 IF (skala_atom_grid) THEN
523 gapw_density_partition = native_skala_gapw_density_partition(my_xc_section)
524 END IF
525 use_virial = ASSOCIATED(virial)
526 IF (use_virial) use_virial = my_calculate_forces .AND. &
527 virial%pv_calculate .AND. (.NOT. virial%pv_numer)
528
529 IF (myfun == xc_none) THEN
530 exc1 = 0.0_dp
531 my_rho_atom_set(:)%exc_h = 0.0_dp
532 my_rho_atom_set(:)%exc_s = 0.0_dp
533 ELSE
534 CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
535 r_val=density_cut)
536 CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
537 r_val=gradient_cut)
538 CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
539 r_val=tau_cut)
540
541 lsd = dft_control%lsd
542 nspins = dft_control%nspins
543 needs = xc_functionals_get_needs(xc_fun_section, &
544 lsd=lsd, &
545 calc_potential=.true.)
546
547 gradient_f = (needs%drho .OR. needs%drho_spin) .OR. skala_atom_grid
548 tau_f = (needs%tau .OR. needs%tau_spin) .OR. skala_atom_grid
549
550 IF (atom_composite_active) THEN
551 IF (lsd) THEN
552 needs%rho_spin = .true.
553 needs%drho_spin = .true.
554 needs%tau_spin = .true.
555 ELSE
556 needs%rho = .true.
557 needs%drho = .true.
558 needs%tau = .true.
559 END IF
560
561 ALLOCATE (composite_atomic_grid_sizes(SIZE(particle_set)), &
562 composite_atom_kind(SIZE(particle_set)), &
563 composite_atom_kind_index(SIZE(particle_set)), &
564 composite_atom_start(SIZE(particle_set)), &
565 composite_atom_end(SIZE(particle_set)), &
566 composite_atom_coords(3, SIZE(particle_set)), &
567 composite_partition_weights(SIZE(particle_set)), &
568 composite_partition_atom_coords(3, SIZE(particle_set)), &
569 composite_distances(SIZE(particle_set)))
570 composite_atomic_grid_sizes = 0_int_8
571 composite_atom_kind = 0
572 composite_atom_kind_index = 0
573 composite_atom_start = 0
574 composite_atom_end = 0
575 DO iatom = 1, SIZE(particle_set)
576 composite_atom_coords(:, iatom) = particle_set(iatom)%r
577 END DO
578 DO ikind = 1, SIZE(atomic_kind_set)
579 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
580 NULLIFY (gth_potential, sgp_potential)
581 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
582 gth_potential=gth_potential, grid_atom=grid_atom, &
583 sgp_potential=sgp_potential, zatom=zatom, zeff=zeff)
584 DO iat = 1, natom
585 iatom = atom_list(iat)
586 composite_atomic_grid_sizes(iatom) = int(grid_atom%nr*grid_atom%ng_sphere, kind=int_8)
587 composite_atom_kind(iatom) = ikind
588 composite_atom_kind_index(iatom) = iat
589 END DO
590 END DO
591 IF (any(composite_atomic_grid_sizes <= 0_int_8)) THEN
592 CALL cp_abort(__location__, &
593 "The atom-composite diagnostic requires a GAPW one-center grid for every atom.")
594 END IF
595
596 composite_local_natom = 0
597 DO ikind = 1, SIZE(atomic_kind_set)
598 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
599 NULLIFY (gth_potential, sgp_potential)
600 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
601 gth_potential=gth_potential, sgp_potential=sgp_potential, &
602 zatom=zatom, zeff=zeff)
603 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
604 composite_local_natom = composite_local_natom + max(0, bo(2) - bo(1) + 1)
605 END DO
606 ALLOCATE (composite_local_atoms(composite_local_natom), &
607 composite_local_grid_sizes(composite_local_natom), &
608 composite_local_atom_coords(3, composite_local_natom))
609 composite_local_atom = 0
610 composite_nflat = 0
611 DO ikind = 1, SIZE(atomic_kind_set)
612 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
613 NULLIFY (gth_potential, sgp_potential)
614 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
615 gth_potential=gth_potential, sgp_potential=sgp_potential, &
616 zatom=zatom, zeff=zeff)
617 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
618 DO iat = bo(1), bo(2)
619 iatom = atom_list(iat)
620 composite_local_atom = composite_local_atom + 1
621 composite_local_atoms(composite_local_atom) = iatom
622 composite_local_grid_sizes(composite_local_atom) = &
623 composite_atomic_grid_sizes(iatom)
624 composite_local_atom_coords(:, composite_local_atom) = &
625 composite_atom_coords(:, iatom)
626 composite_atom_start(iatom) = composite_nflat + 1
627 composite_nflat = composite_nflat + int(composite_atomic_grid_sizes(iatom))
628 composite_atom_end(iatom) = composite_nflat
629 END DO
630 END DO
631 cpassert(composite_local_atom == composite_local_natom)
632 ALLOCATE (composite_density(composite_nflat, 2), &
633 composite_grad(composite_nflat, 3, 2), &
634 composite_kin(composite_nflat, 2), &
635 composite_grid_atom(composite_nflat), &
636 composite_smooth_density_cache(composite_nflat, 2), &
637 composite_smooth_gradient_cache(composite_nflat, 3, 2), &
638 composite_smooth_kin_cache(composite_nflat, 2), &
639 composite_grid_coords(3, composite_nflat), &
640 composite_grid_weights(composite_nflat), &
641 composite_base_grid_weights(composite_nflat), &
642 composite_atomic_grid_weights(composite_nflat))
643 composite_density = 0.0_dp
644 composite_grad = 0.0_dp
645 composite_kin = 0.0_dp
646 DO composite_local_atom = 1, composite_local_natom
647 iatom = composite_local_atoms(composite_local_atom)
648 composite_grid_atom(composite_atom_start(iatom):composite_atom_end(iatom)) = iatom
649 END DO
650 composite_smooth_density_cache = 0.0_dp
651 composite_smooth_gradient_cache = 0.0_dp
652 composite_smooth_kin_cache = 0.0_dp
653 composite_grid_coords = 0.0_dp
654 composite_grid_weights = 0.0_dp
655 composite_base_grid_weights = 0.0_dp
656 composite_atomic_grid_weights = 0.0_dp
657
658 CALL qs_rho_get(rho_struct, rho_r=smooth_rho_r, rho_g=smooth_rho_g, &
659 tau_r=smooth_tau_r, rho_g_valid=rho_g_valid, &
660 tau_r_valid=tau_r_valid)
661 cpassert(rho_g_valid)
662 cpassert(tau_r_valid)
663 cpassert(ASSOCIATED(smooth_rho_r))
664 cpassert(ASSOCIATED(smooth_rho_g))
665 cpassert(ASSOCIATED(smooth_tau_r))
666 CALL prepare_native_grid_cache(qs_env%native_grid_cache, smooth_rho_r(1)%pw_grid, &
667 cell, composite_nflat, image_partition_atom_composite)
668 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
669 CALL section_vals_val_get(my_xc_section, "XC_GRID%XC_DERIV", &
670 i_val=xc_deriv_method_id)
671 CALL section_vals_val_get(my_xc_section, "XC_GRID%XC_SMOOTH_RHO", &
672 i_val=xc_rho_smooth_id)
673 CALL xc_rho_set_create(smooth_rho_set, smooth_rho_r(1)%pw_grid%bounds_local, &
674 rho_cutoff=section_get_rval(my_xc_section, "density_cutoff"), &
675 drho_cutoff=section_get_rval(my_xc_section, "gradient_cutoff"), &
676 tau_cutoff=section_get_rval(my_xc_section, "tau_cutoff"))
677 CALL xc_rho_set_update(smooth_rho_set, smooth_rho_r, smooth_rho_g, smooth_tau_r, needs, &
678 xc_deriv_method_id, xc_rho_smooth_id, auxbas_pw_pool)
679 IF (lsd) THEN
680 CALL xc_rho_set_get(smooth_rho_set, rhoa=smooth_rhoa, rhob=smooth_rhob, &
681 drhoa=smooth_drhoa, drhob=smooth_drhob, &
682 tau_a=smooth_tau_a, tau_b=smooth_tau_b)
683 CALL gather_native_grid_field(smooth_rhoa, smooth_rho_r(1)%pw_grid, para_env, &
684 composite_smooth_rhoa)
685 CALL gather_native_grid_field(smooth_rhob, smooth_rho_r(1)%pw_grid, para_env, &
686 composite_smooth_rhob)
687 CALL gather_native_grid_field(smooth_tau_a, smooth_rho_r(1)%pw_grid, para_env, &
688 composite_smooth_tau_a)
689 CALL gather_native_grid_field(smooth_tau_b, smooth_rho_r(1)%pw_grid, para_env, &
690 composite_smooth_tau_b)
691 DO idir = 1, 3
692 CALL gather_native_grid_field(smooth_drhoa(idir)%array, &
693 smooth_rho_r(1)%pw_grid, para_env, &
694 composite_smooth_drhoa(idir)%array)
695 CALL gather_native_grid_field(smooth_drhob(idir)%array, &
696 smooth_rho_r(1)%pw_grid, para_env, &
697 composite_smooth_drhob(idir)%array)
698 END DO
699 ELSE
700 CALL xc_rho_set_get(smooth_rho_set, rho=smooth_rho, drho=smooth_drho, &
701 tau=smooth_tau)
702 CALL gather_native_grid_field(smooth_rho, smooth_rho_r(1)%pw_grid, para_env, &
703 composite_smooth_rho)
704 CALL gather_native_grid_field(smooth_tau, smooth_rho_r(1)%pw_grid, para_env, &
705 composite_smooth_tau)
706 DO idir = 1, 3
707 CALL gather_native_grid_field(smooth_drho(idir)%array, &
708 smooth_rho_r(1)%pw_grid, para_env, &
709 composite_smooth_drho(idir)%array)
710 END DO
711 END IF
712 END IF
713
714 ! Initialize energy contribution from the one center XC terms to zero
715 exc1 = 0.0_dp
716
717 ! Nullify some pointers for work-arrays
718 NULLIFY (rho_h, drho_h, rho_s, drho_s, weight_h, weight_s)
719 NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
720 NULLIFY (tau_h, tau_s)
721 NULLIFY (vtau_h, vtau_s)
722
723 ! Here starts the loop over all the atoms
724
725 DO ikind = 1, SIZE(atomic_kind_set)
726 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
727 NULLIFY (gth_potential, sgp_potential)
728 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
729 gth_potential=gth_potential, harmonics=harmonics, &
730 grid_atom=grid_atom, sgp_potential=sgp_potential, &
731 zatom=zatom, zeff=zeff)
732 one_center_kind = .NOT. direct_valence_atom_composite
733 IF (one_center_kind) THEN
734 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
735 one_center_kind = paw_atom
736 IF (skala_atom_grid) THEN
737 one_center_kind = native_skala_uses_one_center_kind( &
738 paw_atom, gapw_representation, &
739 ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential), &
740 zeff, zatom)
741 END IF
742 END IF
743 IF (.NOT. one_center_kind .AND. .NOT. atom_composite_active) cycle
744
745 nr = grid_atom%nr
746 na = grid_atom%ng_sphere
747
748 IF (one_center_kind) THEN
749 ! Prepare the structures needed to calculate and store the one-center XC derivatives.
750
751 ! Array dimension: here anly one dimensional arrays are used,
752 ! i.e. only the first column of deriv_data is read.
753 ! The other to dimensions are set to size equal 1
754 bounds(1:2, 1:3) = 1
755 bounds(2, 1) = na
756 bounds(2, 2) = nr
757
758 ! set integration weights
759 weight_h => grid_atom%weight
760 weight_s => grid_atom%weight
761 IF (accint) THEN
762 on = dft_control%qs_control%gapw_control%oweights
763 alpha = dft_control%qs_control%gapw_control%aw(ikind)
764 IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
765 IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
766 END IF
767 IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
768 ALLOCATE (grid_atom%gapw_weight_s(na, nr))
769 ALLOCATE (fw(nr))
770 CALL calc_weight_function(fw, grid_atom%rad2, on, alpha)
771 DO ir = 1, nr
772 agr = 1.0_dp - fw(ir)
773 grid_atom%gapw_weight_s(:, ir) = agr*grid_atom%weight(:, ir)
774 END DO
775 DEALLOCATE (fw)
776 grid_atom%gapw_weight_alpha = alpha
777 END IF
778 weight_s => grid_atom%gapw_weight_s
779 END IF
780
781 ! create a place where to put the derivatives
782 CALL xc_dset_create(deriv_set, local_bounds=bounds)
783 ! create the place where to store the argument for the functionals
784 CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
785 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
786 CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
787 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
788
789 ! allocate the required 3d arrays where to store rho and drho
790 CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
791 CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
792
793 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
794 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
795 CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
796 CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
797 !
798 IF (gradient_f) THEN
799 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
800 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
801 CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
802 CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
803 END IF
804
805 IF (tau_f) THEN
806 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
807 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
808 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
809 CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
810 CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
811 END IF
812
813 ! NLCC for separate hard and soft one-center densities.
814 donlcc = .false.
815 IF (nlcc) THEN
816 NULLIFY (rho_nlcc)
817 rho_nlcc => my_kind_set(ikind)%nlcc_pot
818 IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
819 END IF
820 END IF
821
822 ! Distribute the atoms of this kind
823
824 num_pe = para_env%num_pe
825 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
826
827 DO iat = bo(1), bo(2)
828 iatom = atom_list(iat)
829
830 IF (one_center_kind) THEN
831 my_rho_atom_set(iatom)%exc_h = 0.0_dp
832 my_rho_atom_set(iatom)%exc_s = 0.0_dp
833
834 rho_atom => my_rho_atom_set(iatom)
835 rho_h = 0.0_dp
836 rho_s = 0.0_dp
837 IF (gradient_f) THEN
838 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
839 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
840 rho_rad_s=r_s, drho_rad_h=dr_h, &
841 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
842 rho_rad_s_d=r_s_d)
843 drho_h = 0.0_dp
844 drho_s = 0.0_dp
845 ELSE
846 NULLIFY (r_h, r_s)
847 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
848 rho_d = 0.0_dp
849 END IF
850 IF (tau_f) THEN
851 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
852 ELSE
853 tau_d = 0.0_dp
854 END IF
855
856 DO ir = 1, nr
857 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
858 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
859 r_h_d, r_s_d, drho_h, drho_s)
860 IF (donlcc) THEN
861 CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
862 ir, rho_nlcc(:, 1), rho_h, rho_s, &
863 rho_nlcc(:, 2), drho_h, drho_s)
864 END IF
865 END DO
866 END IF
867
868 IF (atom_composite_active) THEN
869 IF (image_partition_atom_composite) THEN
870 CALL build_periodic_atom_image_layout( &
871 composite_atom_coords, cell, iatom, composite_image_periodicity, &
872 composite_partition_image_coords, composite_partition_target_image)
873 CALL build_periodic_atom_image_layout( &
874 composite_atom_coords(:, iatom:iatom), cell, 1, &
875 composite_image_periodicity, composite_descriptor_image_coords, &
876 composite_descriptor_target_image)
877 END IF
878!$OMP PARALLEL DO COLLAPSE(2) IF (image_partition_atom_composite) SCHEDULE(STATIC) DEFAULT(NONE) &
879!$OMP PRIVATE(composite_row, composite_point, composite_smooth_density_value, &
880!$OMP composite_smooth_gradient_value, composite_smooth_kin_value, &
881!$OMP descriptor_window_weight, idir, ispin, &
882!$OMP gth_potential, nlcc_density, nlcc_gradient, nlcc_hessian, nlcc_spin_factor, &
883!$OMP interpolation_stencil, partition_scale, partition_weight, sgp_potential, source_atom) &
884!$OMP SHARED(atom_composite_reference, cell, composite_atom_coords, composite_atom_kind, &
885!$OMP composite_atom_start, composite_atomic_grid_weights, composite_base_grid_weights, &
886!$OMP composite_density, composite_descriptor_image_coords, &
887!$OMP composite_descriptor_target_image, composite_distances, composite_grad, &
888!$OMP composite_grid_coords, composite_grid_weights, composite_kin, &
889!$OMP composite_partition_atom_coords, composite_partition_image_coords, &
890!$OMP composite_partition_target_image, composite_partition_weights, &
891!$OMP composite_smooth_density_cache, composite_smooth_drho, composite_smooth_drhoa, &
892!$OMP composite_smooth_drhob, composite_smooth_gradient_cache, &
893!$OMP composite_smooth_kin_cache, composite_smooth_rho, composite_smooth_rhoa, &
894!$OMP composite_smooth_rhob, composite_smooth_tau, composite_smooth_tau_a, &
895!$OMP composite_smooth_tau_b, drho_h, drho_s, grid_atom, iatom, &
896!$OMP image_partition_atom_composite, lsd, my_kind_set, na, nlcc, nr, one_center_kind, &
897!$OMP particle_set, qs_env, rho_h, rho_s, smooth_rho_r, tau_h, tau_s, &
898!$OMP use_atom_composite_density, use_atom_composite_gradient, use_atom_composite_tau)
899 DO ir = 1, nr
900 DO ia = 1, na
901 composite_row = composite_atom_start(iatom) + (ir - 1)*na + ia - 1
902 composite_point(1) = particle_set(iatom)%r(1) + grid_atom%rad(ir)* &
903 grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
904 composite_point(2) = particle_set(iatom)%r(2) + grid_atom%rad(ir)* &
905 grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
906 composite_point(3) = particle_set(iatom)%r(3) + &
907 grid_atom%rad(ir)*grid_atom%cos_pol(ia)
908 composite_grid_coords(:, composite_row) = composite_point
909 IF (image_partition_atom_composite) THEN
910 CALL periodic_atom_image_partition_from_layout( &
911 composite_point, composite_partition_image_coords, &
912 composite_partition_target_image, partition_weight)
913 ! The self-image partition defines a smooth atom-centered periodic
914 ! descriptor domain without truncating it at neighboring atoms.
915 CALL periodic_atom_image_partition_from_layout( &
916 composite_point, composite_descriptor_image_coords, &
917 composite_descriptor_target_image, descriptor_window_weight)
918 partition_scale = smooth_partition_atomic_weight_scale( &
919 descriptor_window_weight)
920 ELSE
921 CALL smooth_atom_partition( &
922 composite_point, composite_atom_coords, cell, &
923 composite_partition_weights, composite_partition_atom_coords, &
924 composite_distances)
925 partition_weight = composite_partition_weights(iatom)
926 partition_scale = 1.0_dp
927 END IF
928 composite_base_grid_weights(composite_row) = grid_atom%weight(ia, ir)
929 composite_atomic_grid_weights(composite_row) = &
930 composite_base_grid_weights(composite_row)*partition_scale
931 composite_grid_weights(composite_row) = &
932 composite_base_grid_weights(composite_row)* &
933 partition_weight
934 IF (.NOT. fetch_native_grid_stencil(qs_env%native_grid_cache, composite_row, &
935 composite_point, interpolation_stencil)) THEN
936 CALL create_native_grid_interpolation_stencil( &
937 interpolation_stencil, smooth_rho_r(1)%pw_grid, cell, composite_point, &
938 image_partition_atom_composite)
939 CALL store_native_grid_stencil(qs_env%native_grid_cache, composite_row, &
940 composite_point, interpolation_stencil)
941 END IF
942 IF (lsd) THEN
943 CALL interpolate_native_grid_fields( &
944 composite_smooth_rhoa, composite_smooth_drhoa(1)%array, &
945 composite_smooth_drhoa(2)%array, composite_smooth_drhoa(3)%array, &
946 composite_smooth_tau_a, interpolation_stencil, &
947 composite_smooth_density_value(1), &
948 composite_smooth_gradient_value(:, 1), &
949 composite_smooth_kin_value(1))
950 CALL interpolate_native_grid_fields( &
951 composite_smooth_rhob, composite_smooth_drhob(1)%array, &
952 composite_smooth_drhob(2)%array, composite_smooth_drhob(3)%array, &
953 composite_smooth_tau_b, interpolation_stencil, &
954 composite_smooth_density_value(2), &
955 composite_smooth_gradient_value(:, 2), &
956 composite_smooth_kin_value(2))
957 composite_smooth_density_cache(composite_row, :) = &
958 composite_smooth_density_value
959 composite_smooth_gradient_cache(composite_row, :, :) = &
960 composite_smooth_gradient_value
961 composite_smooth_kin_cache(composite_row, :) = &
962 composite_smooth_kin_value
963 DO ispin = 1, 2
964 composite_density(composite_row, ispin) = &
965 composite_smooth_density_value(ispin)
966 composite_grad(composite_row, :, ispin) = &
967 composite_smooth_gradient_value(:, ispin)
968 composite_kin(composite_row, ispin) = &
969 composite_smooth_kin_value(ispin)
970 IF (one_center_kind .AND. use_atom_composite_density) THEN
971 composite_density(composite_row, ispin) = &
972 composite_density(composite_row, ispin) + &
973 rho_h(ia, ir, ispin) - rho_s(ia, ir, ispin)
974 END IF
975 IF (one_center_kind .AND. use_atom_composite_gradient) THEN
976 DO idir = 1, 3
977 composite_grad(composite_row, idir, ispin) = &
978 composite_grad(composite_row, idir, ispin) + &
979 drho_h(idir, ia, ir, ispin) - drho_s(idir, ia, ir, ispin)
980 END DO
981 END IF
982 IF (one_center_kind .AND. use_atom_composite_tau) THEN
983 composite_kin(composite_row, ispin) = &
984 composite_kin(composite_row, ispin) + &
985 tau_h(ia, ir, ispin) - tau_s(ia, ir, ispin)
986 END IF
987 END DO
988 ELSE
989 CALL interpolate_native_grid_fields( &
990 composite_smooth_rho, composite_smooth_drho(1)%array, &
991 composite_smooth_drho(2)%array, composite_smooth_drho(3)%array, &
992 composite_smooth_tau, interpolation_stencil, &
993 composite_smooth_density_value(1), &
994 composite_smooth_gradient_value(:, 1), &
995 composite_smooth_kin_value(1))
996 composite_smooth_density_cache(composite_row, 1) = &
997 composite_smooth_density_value(1)
998 composite_smooth_gradient_cache(composite_row, :, 1) = &
999 composite_smooth_gradient_value(:, 1)
1000 composite_smooth_kin_cache(composite_row, 1) = &
1001 composite_smooth_kin_value(1)
1002 composite_density(composite_row, :) = &
1003 0.5_dp*composite_smooth_density_value(1)
1004 DO idir = 1, 3
1005 composite_grad(composite_row, idir, :) = &
1006 0.5_dp*composite_smooth_gradient_value(idir, 1)
1007 END DO
1008 composite_kin(composite_row, :) = &
1009 0.5_dp*composite_smooth_kin_value(1)
1010 IF (one_center_kind .AND. use_atom_composite_density) THEN
1011 composite_density(composite_row, :) = &
1012 composite_density(composite_row, :) + &
1013 0.5_dp*(rho_h(ia, ir, 1) - rho_s(ia, ir, 1))
1014 END IF
1015 IF (one_center_kind .AND. use_atom_composite_gradient) THEN
1016 DO idir = 1, 3
1017 composite_grad(composite_row, idir, :) = &
1018 composite_grad(composite_row, idir, :) + &
1019 0.5_dp*(drho_h(idir, ia, ir, 1) - &
1020 drho_s(idir, ia, ir, 1))
1021 END DO
1022 END IF
1023 IF (one_center_kind .AND. use_atom_composite_tau) THEN
1024 composite_kin(composite_row, :) = composite_kin(composite_row, :) + &
1025 0.5_dp*(tau_h(ia, ir, 1) - &
1026 tau_s(ia, ir, 1))
1027 END IF
1028 END IF
1029 IF (atom_composite_reference .AND. nlcc) THEN
1030 nlcc_spin_factor = merge(1.0_dp, 0.5_dp, lsd)
1031 DO source_atom = 1, SIZE(particle_set)
1032 NULLIFY (gth_potential, sgp_potential)
1033 CALL get_qs_kind(my_kind_set(composite_atom_kind(source_atom)), &
1034 gth_potential=gth_potential, &
1035 sgp_potential=sgp_potential)
1036 CALL evaluate_nlcc_primitive_fields( &
1037 composite_point, particle_set(source_atom)%r, &
1038 gth_potential, sgp_potential, nlcc_density, &
1039 nlcc_gradient, nlcc_hessian)
1040 composite_density(composite_row, :) = &
1041 composite_density(composite_row, :) + &
1042 nlcc_spin_factor*nlcc_density
1043 DO idir = 1, 3
1044 composite_grad(composite_row, idir, :) = &
1045 composite_grad(composite_row, idir, :) + &
1046 nlcc_spin_factor*nlcc_gradient(idir)
1047 END DO
1048 END DO
1049 END IF
1050 END DO
1051 END DO
1052!$OMP END PARALLEL DO
1053 IF (image_partition_atom_composite) THEN
1054 DEALLOCATE (composite_descriptor_image_coords, composite_partition_image_coords)
1055 END IF
1056 cpassert(nr*na == composite_atom_end(iatom) - composite_atom_start(iatom) + 1)
1057 END IF
1058
1059 IF (.NOT. one_center_kind) cycle
1060
1061 DO ir = 1, nr
1062 IF (tau_f) THEN
1063 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
1064 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
1065 ELSE IF (gradient_f) THEN
1066 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
1067 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
1068 ELSE
1069 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
1070 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
1071 END IF
1072 END DO
1073
1074 evaluate_hard = .true.
1075 evaluate_soft = .true.
1076 skala_atom_force_h = 0.0_dp
1077 skala_atom_force_s = 0.0_dp
1078 skala_atom_virial_h = 0.0_dp
1079 skala_atom_virial_s = 0.0_dp
1080 IF (skala_atom_grid) THEN
1081 SELECT CASE (gapw_density_partition)
1082 CASE (skala_gapw_density_partition_hard_minus_soft)
1083 CONTINUE
1084 CASE (skala_gapw_density_partition_hard_only)
1085 evaluate_soft = .false.
1086 CASE (skala_gapw_density_partition_soft_only)
1087 evaluate_hard = .false.
1088 CASE (skala_gapw_density_partition_none)
1089 evaluate_hard = .false.
1090 evaluate_soft = .false.
1091 CASE DEFAULT
1092 CALL cp_abort(__location__, &
1093 "Unknown GAUXC%NATIVE_GRID_GAPW_DENSITY_PARTITION value.")
1094 END SELECT
1095 END IF
1096 IF (atom_composite_reference) THEN
1097 evaluate_hard = .false.
1098 evaluate_soft = .false.
1099 END IF
1100
1101 !-------------------!
1102 ! hard atom density !
1103 !-------------------!
1104 CALL xc_dset_zero_all(deriv_set)
1105 IF (.NOT. evaluate_hard) THEN
1106 exc_h = 0.0_dp
1107 IF (.NOT. energy_only) THEN
1108 vxc_h = 0.0_dp
1109 IF (ASSOCIATED(vxg_h)) vxg_h = 0.0_dp
1110 IF (ASSOCIATED(vtau_h)) vtau_h = 0.0_dp
1111 END IF
1112 ELSE IF (skala_atom_grid) THEN
1113 CALL skala_gapw_atom_vxc_of_r( &
1114 my_xc_section, grid_atom, para_env, particle_set(iatom)%r, &
1115 rho_h, drho_h, tau_h, weight_h, lsd, nspins, na, nr, &
1116 exc_h, vxc_h, vxg_h, vtau_h, energy_only=energy_only, &
1117 atom_force=skala_atom_force_h, atom_virial=skala_atom_virial_h)
1118 ELSE
1119 CALL vxc_of_r_new(xc_fun_section, rho_set_h, deriv_set, 1, needs, weight_h, &
1120 lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h, energy_only=energy_only, &
1121 adiabatic_rescale_factor=my_adiabatic_rescale_factor)
1122 END IF
1123 rho_atom%exc_h = rho_atom%exc_h + exc_h
1124
1125 !-------------------!
1126 ! soft atom density !
1127 !-------------------!
1128 CALL xc_dset_zero_all(deriv_set)
1129 IF (.NOT. evaluate_soft) THEN
1130 exc_s = 0.0_dp
1131 IF (.NOT. energy_only) THEN
1132 vxc_s = 0.0_dp
1133 IF (ASSOCIATED(vxg_s)) vxg_s = 0.0_dp
1134 IF (ASSOCIATED(vtau_s)) vtau_s = 0.0_dp
1135 END IF
1136 ELSE IF (skala_atom_grid) THEN
1137 CALL skala_gapw_atom_vxc_of_r( &
1138 my_xc_section, grid_atom, para_env, particle_set(iatom)%r, &
1139 rho_s, drho_s, tau_s, weight_s, lsd, nspins, na, nr, &
1140 exc_s, vxc_s, vxg_s, vtau_s, energy_only=energy_only, &
1141 atom_force=skala_atom_force_s, atom_virial=skala_atom_virial_s)
1142 ELSE
1143 CALL vxc_of_r_new(xc_fun_section, rho_set_s, deriv_set, 1, needs, weight_s, &
1144 lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s, energy_only=energy_only, &
1145 adiabatic_rescale_factor=my_adiabatic_rescale_factor)
1146 END IF
1147 rho_atom%exc_s = rho_atom%exc_s + exc_s
1148
1149 ! Add contributions to the exc energy
1150
1151 exc1 = exc1 + rho_atom%exc_h - rho_atom%exc_s
1152 IF (skala_atom_grid .AND. my_calculate_forces .AND. ASSOCIATED(force)) THEN
1153 force(ikind)%rho_elec(:, iat) = force(ikind)%rho_elec(:, iat) + &
1154 skala_atom_force_h - skala_atom_force_s
1155 END IF
1156 IF (skala_atom_grid .AND. use_virial) THEN
1157 skala_atom_virial = skala_atom_virial_h - skala_atom_virial_s
1158 DO idir = 1, 3
1159 DO jdir = 1, 3
1160 virial%pv_gapw(idir, jdir) = virial%pv_gapw(idir, jdir) + &
1161 skala_atom_virial(idir, jdir)
1162 virial%pv_virial(idir, jdir) = virial%pv_virial(idir, jdir) + &
1163 skala_atom_virial(idir, jdir)
1164 END DO
1165 END DO
1166 END IF
1167
1168 ! Integration to get the matrix elements relative to the vxc_atom
1169 ! here the products with the primitives is done: gaVxcgb
1170 ! internal transformation to get the integral in cartesian Gaussians
1171
1172 IF (.NOT. energy_only) THEN
1173 NULLIFY (int_hh, int_ss)
1174 CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
1175 IF (gradient_f) THEN
1176 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
1177 grid_atom, basis_1c, harmonics, nspins)
1178 ELSE
1179 CALL gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, &
1180 grid_atom, basis_1c, harmonics, nspins)
1181 END IF
1182 IF (tau_f) THEN
1183 CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, tau_basis_cache, nspins)
1184 END IF
1185 END IF ! energy_only
1186 NULLIFY (r_h, r_s, dr_h, dr_s)
1187 END DO ! iat
1188
1189 IF (one_center_kind) THEN
1190 IF (tau_f) CALL release_tau_basis_cache(tau_basis_cache)
1191
1192 CALL xc_dset_release(deriv_set)
1193 CALL xc_rho_set_release(rho_set_h)
1194 CALL xc_rho_set_release(rho_set_s)
1195 END IF
1196 END DO ! ikind
1197
1198 IF (atom_composite_active) THEN
1199 ALLOCATE (composite_cross_density(composite_nflat, 2), &
1200 composite_cross_grad(composite_nflat, 3, 2), &
1201 composite_cross_kin(composite_nflat, 2))
1202 composite_cross_density = 0.0_dp
1203 composite_cross_grad = 0.0_dp
1204 composite_cross_kin = 0.0_dp
1205 composite_cross_cutoff_max = 0.0_dp
1206 IF (.NOT. direct_valence_atom_composite) THEN
1207 DO ikind = 1, SIZE(atomic_kind_set)
1208 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1209 NULLIFY (gth_potential, sgp_potential)
1210 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
1211 gth_potential=gth_potential, harmonics=harmonics, &
1212 grid_atom=grid_atom, sgp_potential=sgp_potential, &
1213 zatom=zatom, zeff=zeff)
1214 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1215 IF (.NOT. native_skala_uses_one_center_kind( &
1216 paw_atom, gapw_representation, &
1217 ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential), &
1218 zeff, zatom)) cycle
1219
1220 CALL replicate_rho_atom_radial( &
1221 para_env, my_rho_atom_set, my_kind_set(ikind), atom_list, natom, nspins)
1222 nr = grid_atom%nr
1223 na = grid_atom%ng_sphere
1224 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
1225 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
1226 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
1227 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
1228 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
1229 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
1230 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
1231
1232 ! The one-center density matrices are already globally reduced. Distribute the
1233 ! overlap work by target atom so that every rank constructs only its model rows.
1234 DO iat = 1, natom
1235 source_atom = atom_list(iat)
1236 rho_atom => my_rho_atom_set(source_atom)
1237 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1238 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s, &
1239 drho_rad_h=dr_h, drho_rad_s=dr_s, &
1240 rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1241 rho_h = 0.0_dp
1242 rho_s = 0.0_dp
1243 drho_h = 0.0_dp
1244 drho_s = 0.0_dp
1245 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
1246 DO ir = 1, nr
1247 CALL calc_rho_angular(grid_atom, harmonics, nspins, .true., &
1248 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
1249 r_h_d, r_s_d, drho_h, drho_s)
1250 END DO
1251
1252 cross_cutoff = gapw_atom_grid_support_radius( &
1253 grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s)
1254 IF (cross_cutoff <= 0.0_dp) cycle
1255 composite_cross_cutoff_max = max(composite_cross_cutoff_max, cross_cutoff)
1256 image_shell = 0
1257 DO idir = 1, 3
1258 IF (cell%perd(idir) == 1) THEN
1259 image_shell(idir) = ceiling( &
1260 cross_cutoff*sqrt(sum(cell%h_inv(idir, :)**2))) + 1
1261 END IF
1262 END DO
1263
1264!$OMP PARALLEL DO DEFAULT(NONE) SCHEDULE(STATIC) &
1265!$OMP PRIVATE(image_lower, image_upper, composite_row, cross_density, cross_density_spatial, &
1266!$OMP cross_displacement, cross_grad, cross_grad_spatial, cross_kin, &
1267!$OMP cross_kin_spatial, fractional, idir, image_i1, image_i2, image_i3, jdir, &
1268!$OMP image_shift, image_translation, target_atom) &
1269!$OMP SHARED(cell, composite_cross_density, composite_cross_grad, composite_cross_kin, &
1270!$OMP composite_grid_atom, composite_grid_coords, composite_nflat, cross_cutoff, &
1271!$OMP drho_h, drho_s, grid_atom, harmonics, image_shell, lsd, nspins, &
1272!$OMP particle_set, rho_h, rho_s, source_atom, tau_h, tau_s, &
1273!$OMP use_atom_composite_density, use_atom_composite_gradient, use_atom_composite_tau)
1274 DO composite_row = 1, composite_nflat
1275 target_atom = composite_grid_atom(composite_row)
1276 fractional = 0.0_dp
1277 DO idir = 1, 3
1278 DO jdir = 1, 3
1279 fractional(idir) = fractional(idir) + cell%h_inv(idir, jdir)* &
1280 (composite_grid_coords(jdir, composite_row) - &
1281 particle_set(source_atom)%r(jdir))
1282 END DO
1283 END DO
1284 CALL atom_grid_image_bounds(cell, fractional, cross_cutoff, image_shell, &
1285 image_lower, image_upper)
1286 DO image_i3 = image_lower(3), image_upper(3)
1287 DO image_i2 = image_lower(2), image_upper(2)
1288 DO image_i1 = image_lower(1), image_upper(1)
1289 image_shift = [image_i1, image_i2, image_i3]
1290 IF (target_atom == source_atom .AND. &
1291 all(image_shift == 0)) cycle
1292 image_translation = matmul( &
1293 cell%hmat, real(image_shift, dp))
1294 cross_displacement = composite_grid_coords(:, composite_row) - &
1295 particle_set(source_atom)%r - &
1296 image_translation
1297 IF (sqrt(sum(cross_displacement**2)) > cross_cutoff) cycle
1298 CALL interpolate_gapw_atom_grid_fields( &
1299 grid_atom, harmonics, cross_displacement, cross_cutoff, nspins, &
1300 rho_h, rho_s, drho_h, drho_s, tau_h, tau_s, &
1301 cross_density, cross_grad, cross_kin, cross_density_spatial, &
1302 cross_grad_spatial, cross_kin_spatial, calculate_spatial=.false.)
1303 IF (lsd) THEN
1304 IF (use_atom_composite_density) THEN
1305 composite_cross_density(composite_row, 1:2) = &
1306 composite_cross_density(composite_row, 1:2) + &
1307 cross_density(1:2)
1308 END IF
1309 IF (use_atom_composite_gradient) THEN
1310 composite_cross_grad(composite_row, :, 1:2) = &
1311 composite_cross_grad(composite_row, :, 1:2) + &
1312 cross_grad(:, 1:2)
1313 END IF
1314 IF (use_atom_composite_tau) THEN
1315 composite_cross_kin(composite_row, 1:2) = &
1316 composite_cross_kin(composite_row, 1:2) + cross_kin(1:2)
1317 END IF
1318 ELSE
1319 IF (use_atom_composite_density) THEN
1320 composite_cross_density(composite_row, :) = &
1321 composite_cross_density(composite_row, :) + &
1322 0.5_dp*cross_density(1)
1323 END IF
1324 IF (use_atom_composite_gradient) THEN
1325 DO idir = 1, 3
1326 composite_cross_grad(composite_row, idir, :) = &
1327 composite_cross_grad(composite_row, idir, :) + &
1328 0.5_dp*cross_grad(idir, 1)
1329 END DO
1330 END IF
1331 IF (use_atom_composite_tau) THEN
1332 composite_cross_kin(composite_row, :) = &
1333 composite_cross_kin(composite_row, :) + &
1334 0.5_dp*cross_kin(1)
1335 END IF
1336 END IF
1337 END DO
1338 END DO
1339 END DO
1340 END DO
1341!$OMP END PARALLEL DO
1342 END DO
1343
1344 CALL release_tau_basis_cache(tau_basis_cache)
1345 END DO
1346 END IF
1347
1348 IF (native_grid_diagnostics) THEN
1349 composite_cross_density_max = maxval(abs(composite_cross_density))
1350 composite_cross_grad_max = maxval(abs(composite_cross_grad))
1351 composite_cross_kin_max = maxval(abs(composite_cross_kin))
1352 CALL para_env%max(composite_cross_cutoff_max)
1353 CALL para_env%max(composite_cross_density_max)
1354 CALL para_env%max(composite_cross_grad_max)
1355 CALL para_env%max(composite_cross_kin_max)
1356 iw = cp_logger_get_default_io_unit()
1357 IF (iw > 0) THEN
1358 WRITE (unit=iw, fmt="(T2,A,4(1X,ES20.12))") &
1359 "SKALA_GPW| Atom-composite cross support/maxima", &
1360 composite_cross_cutoff_max, composite_cross_density_max, &
1361 composite_cross_grad_max, composite_cross_kin_max
1362 END IF
1363 END IF
1364 composite_density(:, :) = composite_density(:, :) + composite_cross_density(:, :)
1365 composite_grad(:, :, :) = composite_grad(:, :, :) + composite_cross_grad(:, :, :)
1366 composite_kin(:, :) = composite_kin(:, :) + composite_cross_kin(:, :)
1367 DEALLOCATE (composite_cross_density, composite_cross_grad, composite_cross_kin)
1368
1369 cpassert(all(composite_grid_weights >= 0.0_dp))
1370 atom_composite_nelec = sum(composite_grid_weights* &
1371 (composite_density(:, 1) + composite_density(:, 2)))
1372 CALL para_env%sum(atom_composite_nelec)
1373 IF (atom_composite_reference .AND. my_calculate_forces) THEN
1374 CALL skala_gapw_atom_composite_energy( &
1375 my_xc_section, para_env, composite_density, composite_grad, composite_kin, &
1376 composite_grid_coords, composite_grid_weights, composite_atomic_grid_weights, &
1377 composite_local_grid_sizes, composite_local_atom_coords, atom_composite_exc, &
1378 composite_density_grad, composite_grad_grad, composite_kin_grad, &
1379 composite_grid_coord_grad, composite_grid_weight_grad, &
1380 composite_atomic_grid_weight_grad, composite_atom_coord_grad)
1381 ALLOCATE (composite_cross_force(3, SIZE(particle_set)), &
1382 composite_explicit_force(3, SIZE(particle_set)), &
1383 composite_grid_coord_force(3, SIZE(particle_set)), &
1384 composite_model_atom_force(3, SIZE(particle_set)), &
1385 composite_moving_smooth_force(3, SIZE(particle_set)), &
1386 composite_nlcc_center_force(3, SIZE(particle_set)), &
1387 composite_nlcc_target_force(3, SIZE(particle_set)), &
1388 composite_partition_force(3, SIZE(particle_set)), &
1389 composite_partition_included(SIZE(particle_set)), &
1390 composite_partition_datom(3, SIZE(particle_set), SIZE(particle_set)), &
1391 composite_partition_dstrain(3, 3, SIZE(particle_set)))
1392 composite_cross_force = 0.0_dp
1393 composite_cross_image_virial = 0.0_dp
1394 composite_model_atom_force = 0.0_dp
1395 composite_grid_coord_force = 0.0_dp
1396 composite_moving_smooth_force = 0.0_dp
1397 composite_nlcc_center_force = 0.0_dp
1398 composite_nlcc_target_force = 0.0_dp
1399 composite_partition_force = 0.0_dp
1400 composite_explicit_virial = 0.0_dp
1401 composite_feature_virial = 0.0_dp
1402 composite_interpolation_virial = 0.0_dp
1403 composite_partition_strain_virial = 0.0_dp
1404!$OMP PARALLEL DO IF (image_partition_atom_composite) SCHEDULE(STATIC) DEFAULT(NONE) &
1405!$OMP PRIVATE(composite_local_atom, composite_row, composite_smooth_gradient_value, &
1406!$OMP iatom, idir, ispin, jdir) &
1407!$OMP SHARED(cell, composite_atom_end, composite_atom_start, composite_grad_grad, &
1408!$OMP composite_grid_coords, composite_local_atoms, composite_local_natom, &
1409!$OMP composite_smooth_drho, composite_smooth_drhoa, composite_smooth_drhob, &
1410!$OMP image_partition_atom_composite, lsd, smooth_rho_r) &
1411!$OMP REDUCTION(+:composite_feature_virial)
1412 DO composite_local_atom = 1, composite_local_natom
1413 iatom = composite_local_atoms(composite_local_atom)
1414 DO composite_row = composite_atom_start(iatom), composite_atom_end(iatom)
1415 IF (lsd) THEN
1416 DO jdir = 1, 3
1417 composite_smooth_gradient_value(jdir, 1) = &
1418 interpolate_native_grid( &
1419 composite_smooth_drhoa(jdir)%array, smooth_rho_r(1)%pw_grid, &
1420 cell, composite_grid_coords(:, composite_row), &
1421 image_partition_atom_composite)
1422 composite_smooth_gradient_value(jdir, 2) = &
1423 interpolate_native_grid( &
1424 composite_smooth_drhob(jdir)%array, smooth_rho_r(1)%pw_grid, &
1425 cell, composite_grid_coords(:, composite_row), &
1426 image_partition_atom_composite)
1427 END DO
1428 ELSE
1429 DO jdir = 1, 3
1430 composite_smooth_gradient_value(jdir, :) = 0.5_dp* &
1431 interpolate_native_grid( &
1432 composite_smooth_drho(jdir)%array, smooth_rho_r(1)%pw_grid, &
1433 cell, composite_grid_coords(:, composite_row), &
1434 image_partition_atom_composite)
1435 END DO
1436 END IF
1437 DO ispin = 1, 2
1438 DO idir = 1, 3
1439 DO jdir = 1, 3
1440 composite_feature_virial(jdir, idir) = &
1441 composite_feature_virial(jdir, idir) - &
1442 composite_grad_grad(composite_row, idir, ispin)* &
1443 composite_smooth_gradient_value(jdir, ispin)
1444 END DO
1445 END DO
1446 END DO
1447 END DO
1448 END DO
1449!$OMP END PARALLEL DO
1450!$OMP PARALLEL IF (image_partition_atom_composite) DEFAULT(NONE) &
1451!$OMP PRIVATE(composite_local_atom, composite_row, descriptor_window_adjoint, &
1452!$OMP descriptor_window_weight, gth_potential, iatom, idir, ispin, jdir, &
1453!$OMP nlcc_density, nlcc_gradient, nlcc_hessian, nlcc_spatial_derivative, &
1454!$OMP nlcc_spin_factor, partition_adjoint, composite_nlcc_center_force_local, &
1455!$OMP composite_partition_force_local, local_descriptor_datom, &
1456!$OMP local_descriptor_dstrain, local_partition_datom, local_partition_dstrain, &
1457!$OMP local_partition_weight, sgp_potential, source_atom, &
1458!$OMP spatial_derivative, target_atom, target_partition_adjoint) &
1459!$OMP REDUCTION(+:composite_interpolation_virial, composite_partition_strain_virial) &
1460!$OMP SHARED(cell, composite_atom_coord_grad, composite_atom_coords, composite_atom_end, &
1461!$OMP composite_atom_kind, composite_atom_start, composite_atomic_grid_weight_grad, &
1462!$OMP composite_atomic_grid_weights, composite_base_grid_weights, composite_density_grad, &
1463!$OMP composite_grad_grad, composite_grid_coord_force, composite_grid_coord_grad, &
1464!$OMP composite_grid_coords, composite_grid_weight_grad, composite_image_periodicity, &
1465!$OMP composite_kin_grad, composite_local_atoms, composite_local_natom, &
1466!$OMP composite_model_atom_force, composite_moving_smooth_force, &
1467!$OMP composite_nlcc_center_force, composite_nlcc_target_force, composite_partition_datom, &
1468!$OMP composite_partition_dstrain, composite_partition_force, composite_partition_included, &
1469!$OMP composite_partition_weights, composite_smooth_drho, composite_smooth_drhoa, &
1470!$OMP composite_smooth_drhob, composite_smooth_rho, composite_smooth_rhoa, &
1471!$OMP composite_smooth_rhob, composite_smooth_tau, composite_smooth_tau_a, &
1472!$OMP composite_smooth_tau_b, image_partition_atom_composite, lsd, my_kind_set, nlcc, &
1473!$OMP particle_set, smooth_rho_r)
1474 ALLOCATE (composite_nlcc_center_force_local(3, SIZE(particle_set)), &
1475 composite_partition_force_local(3, SIZE(particle_set)), &
1476 local_partition_datom(3, SIZE(particle_set)))
1477 composite_nlcc_center_force_local = 0.0_dp
1478 composite_partition_force_local = 0.0_dp
1479!$OMP DO SCHEDULE(DYNAMIC)
1480 DO composite_local_atom = 1, composite_local_natom
1481 iatom = composite_local_atoms(composite_local_atom)
1482 composite_model_atom_force(:, iatom) = &
1483 composite_atom_coord_grad(:, composite_local_atom)
1484 DO composite_row = composite_atom_start(iatom), composite_atom_end(iatom)
1485 composite_grid_coord_force(:, iatom) = &
1486 composite_grid_coord_force(:, iatom) + &
1487 composite_grid_coord_grad(:, composite_row)
1488 IF (lsd) THEN
1489 spatial_derivative = &
1490 composite_density_grad(composite_row, 1)* &
1491 interpolate_native_grid_gradient( &
1492 composite_smooth_rhoa, smooth_rho_r(1)%pw_grid, cell, &
1493 composite_grid_coords(:, composite_row), &
1494 image_partition_atom_composite) + &
1495 composite_density_grad(composite_row, 2)* &
1496 interpolate_native_grid_gradient( &
1497 composite_smooth_rhob, smooth_rho_r(1)%pw_grid, cell, &
1498 composite_grid_coords(:, composite_row), &
1499 image_partition_atom_composite) + &
1500 composite_kin_grad(composite_row, 1)* &
1501 interpolate_native_grid_gradient( &
1502 composite_smooth_tau_a, smooth_rho_r(1)%pw_grid, cell, &
1503 composite_grid_coords(:, composite_row), &
1504 image_partition_atom_composite) + &
1505 composite_kin_grad(composite_row, 2)* &
1506 interpolate_native_grid_gradient( &
1507 composite_smooth_tau_b, smooth_rho_r(1)%pw_grid, cell, &
1508 composite_grid_coords(:, composite_row), &
1509 image_partition_atom_composite)
1510 DO idir = 1, 3
1511 spatial_derivative = spatial_derivative + &
1512 composite_grad_grad(composite_row, idir, 1)* &
1513 interpolate_native_grid_gradient( &
1514 composite_smooth_drhoa(idir)%array, &
1515 smooth_rho_r(1)%pw_grid, cell, &
1516 composite_grid_coords(:, composite_row), &
1517 image_partition_atom_composite) + &
1518 composite_grad_grad(composite_row, idir, 2)* &
1519 interpolate_native_grid_gradient( &
1520 composite_smooth_drhob(idir)%array, &
1521 smooth_rho_r(1)%pw_grid, cell, &
1522 composite_grid_coords(:, composite_row), &
1523 image_partition_atom_composite)
1524 END DO
1525 ELSE
1526 spatial_derivative = 0.5_dp*sum( &
1527 composite_density_grad(composite_row, :))* &
1528 interpolate_native_grid_gradient( &
1529 composite_smooth_rho, smooth_rho_r(1)%pw_grid, cell, &
1530 composite_grid_coords(:, composite_row), &
1531 image_partition_atom_composite) + &
1532 0.5_dp*sum(composite_kin_grad(composite_row, :))* &
1533 interpolate_native_grid_gradient( &
1534 composite_smooth_tau, smooth_rho_r(1)%pw_grid, cell, &
1535 composite_grid_coords(:, composite_row), &
1536 image_partition_atom_composite)
1537 DO idir = 1, 3
1538 spatial_derivative = spatial_derivative + 0.5_dp*sum( &
1539 composite_grad_grad(composite_row, idir, :))* &
1540 interpolate_native_grid_gradient( &
1541 composite_smooth_drho(idir)%array, &
1542 smooth_rho_r(1)%pw_grid, &
1543 cell, composite_grid_coords(:, composite_row), &
1544 image_partition_atom_composite)
1545 END DO
1546 END IF
1547 DO idir = 1, 3
1548 DO jdir = 1, 3
1549 composite_interpolation_virial(idir, jdir) = &
1550 composite_interpolation_virial(idir, jdir) + &
1551 spatial_derivative(idir)*( &
1552 composite_grid_coords(jdir, composite_row) - &
1553 particle_set(iatom)%r(jdir))
1554 END DO
1555 END DO
1556 composite_moving_smooth_force(:, iatom) = &
1557 composite_moving_smooth_force(:, iatom) + spatial_derivative
1558 IF (nlcc) THEN
1559 nlcc_spin_factor = merge(1.0_dp, 0.5_dp, lsd)
1560 DO source_atom = 1, SIZE(particle_set)
1561 NULLIFY (gth_potential, sgp_potential)
1562 CALL get_qs_kind(my_kind_set(composite_atom_kind(source_atom)), &
1563 gth_potential=gth_potential, &
1564 sgp_potential=sgp_potential)
1565 CALL evaluate_nlcc_primitive_fields( &
1566 composite_grid_coords(:, composite_row), &
1567 particle_set(source_atom)%r, gth_potential, sgp_potential, &
1568 nlcc_density, nlcc_gradient, nlcc_hessian)
1569 nlcc_spatial_derivative = 0.0_dp
1570 DO ispin = 1, 2
1571 nlcc_spatial_derivative = nlcc_spatial_derivative + &
1572 nlcc_spin_factor*composite_density_grad(composite_row, ispin)* &
1573 nlcc_gradient
1574 DO idir = 1, 3
1575 DO jdir = 1, 3
1576 nlcc_spatial_derivative(jdir) = &
1577 nlcc_spatial_derivative(jdir) + nlcc_spin_factor* &
1578 composite_grad_grad(composite_row, idir, ispin)* &
1579 nlcc_hessian(idir, jdir)
1580 END DO
1581 END DO
1582 END DO
1583 composite_nlcc_target_force(:, iatom) = &
1584 composite_nlcc_target_force(:, iatom) + nlcc_spatial_derivative
1585 composite_nlcc_center_force_local(:, source_atom) = &
1586 composite_nlcc_center_force_local(:, source_atom) - nlcc_spatial_derivative
1587 END DO
1588 END IF
1589 IF (image_partition_atom_composite) THEN
1590 CALL periodic_atom_image_partition( &
1591 composite_grid_coords(:, composite_row), composite_atom_coords, cell, &
1592 iatom, local_partition_weight, local_partition_datom, &
1593 local_partition_dstrain, composite_image_periodicity)
1594 CALL periodic_atom_image_partition( &
1595 composite_grid_coords(:, composite_row), &
1596 composite_atom_coords(:, iatom:iatom), cell, 1, &
1597 descriptor_window_weight, local_descriptor_datom, &
1598 local_descriptor_dstrain, composite_image_periodicity)
1599 ! Its atom derivative cancels against the moving target grid; periodic
1600 ! image strain remains an explicit contribution to the virial.
1601 target_partition_adjoint = composite_base_grid_weights(composite_row)* &
1602 composite_grid_weight_grad(composite_row)
1603 descriptor_window_adjoint = composite_base_grid_weights(composite_row)* &
1604 composite_atomic_grid_weight_grad(composite_row)* &
1605 smooth_partition_atomic_weight_scale_derivative( &
1606 descriptor_window_weight)
1607 DO target_atom = 1, SIZE(particle_set)
1608 composite_partition_force_local(:, target_atom) = &
1609 composite_partition_force_local(:, target_atom) + &
1610 target_partition_adjoint*local_partition_datom(:, target_atom)
1611 END DO
1612 composite_partition_force_local(:, iatom) = &
1613 composite_partition_force_local(:, iatom) - target_partition_adjoint* &
1614 sum(local_partition_datom, dim=2)
1615 composite_partition_strain_virial = &
1616 composite_partition_strain_virial - target_partition_adjoint* &
1617 local_partition_dstrain - descriptor_window_adjoint* &
1618 local_descriptor_dstrain
1619 ELSE
1620 CALL skala_gpw_smooth_partition_derivatives( &
1621 composite_grid_coords(:, composite_row), composite_atom_coords, cell, &
1622 composite_partition_weights, composite_partition_included, &
1623 composite_partition_datom, composite_partition_dstrain)
1624 partition_adjoint = composite_grid_weight_grad(composite_row)* &
1625 composite_atomic_grid_weights(composite_row)
1626 DO target_atom = 1, SIZE(particle_set)
1627 composite_partition_force_local(:, target_atom) = &
1628 composite_partition_force_local(:, target_atom) + &
1629 partition_adjoint* &
1630 composite_partition_datom(:, target_atom, iatom)
1631 END DO
1632 composite_partition_force_local(:, iatom) = &
1633 composite_partition_force_local(:, iatom) - partition_adjoint* &
1634 sum(composite_partition_datom(:, :, iatom), dim=2)
1635 END IF
1636 END DO
1637 END DO
1638!$OMP END DO
1639!$OMP CRITICAL(skala_atom_composite_force_reduction)
1640 composite_nlcc_center_force(:, :) = composite_nlcc_center_force(:, :) + &
1641 composite_nlcc_center_force_local
1642 composite_partition_force(:, :) = composite_partition_force(:, :) + &
1643 composite_partition_force_local
1644!$OMP END CRITICAL(skala_atom_composite_force_reduction)
1645 DEALLOCATE (composite_nlcc_center_force_local, composite_partition_force_local, &
1646 local_partition_datom)
1647!$OMP END PARALLEL
1648 ! The target atom grids can overlap augmentation regions of other atoms
1649 ! and periodic images. Differentiate the same discrete interpolation used
1650 ! in the forward composite fields, including the explicit image strain.
1651 IF (.NOT. direct_valence_atom_composite) THEN
1652 DO ikind = 1, SIZE(atomic_kind_set)
1653 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
1654 NULLIFY (gth_potential, sgp_potential)
1655 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
1656 gth_potential=gth_potential, harmonics=harmonics, &
1657 grid_atom=grid_atom, sgp_potential=sgp_potential, &
1658 zatom=zatom, zeff=zeff)
1659 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
1660 IF (.NOT. native_skala_uses_one_center_kind( &
1661 paw_atom, gapw_representation, &
1662 ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential), &
1663 zeff, zatom)) cycle
1664
1665 nr = grid_atom%nr
1666 na = grid_atom%ng_sphere
1667 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
1668 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
1669 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
1670 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
1671 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
1672 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
1673 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
1674
1675 DO iat = 1, natom
1676 source_atom = atom_list(iat)
1677 rho_atom => my_rho_atom_set(source_atom)
1678 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
1679 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s, &
1680 drho_rad_h=dr_h, drho_rad_s=dr_s, &
1681 rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d)
1682 rho_h = 0.0_dp
1683 rho_s = 0.0_dp
1684 drho_h = 0.0_dp
1685 drho_s = 0.0_dp
1686 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
1687 DO ir = 1, nr
1688 CALL calc_rho_angular(grid_atom, harmonics, nspins, .true., &
1689 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
1690 r_h_d, r_s_d, drho_h, drho_s)
1691 END DO
1692
1693 cross_cutoff = gapw_atom_grid_support_radius( &
1694 grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s)
1695 IF (cross_cutoff <= 0.0_dp) cycle
1696 image_shell = 0
1697 DO idir = 1, 3
1698 IF (cell%perd(idir) == 1) THEN
1699 image_shell(idir) = ceiling( &
1700 cross_cutoff*sqrt(sum(cell%h_inv(idir, :)**2))) + 1
1701 END IF
1702 END DO
1703
1704!$OMP PARALLEL DEFAULT(NONE) &
1705!$OMP PRIVATE(image_lower, image_upper, composite_cross_force_local, composite_cross_image_virial_local, &
1706!$OMP composite_row, cross_density, cross_density_adjoint, cross_density_spatial, &
1707!$OMP cross_displacement, cross_grad, cross_grad_adjoint, cross_grad_spatial, &
1708!$OMP cross_kin, cross_kin_adjoint, cross_kin_spatial, cross_spatial_derivative, &
1709!$OMP fractional, idir, image_i1, image_i2, image_i3, image_shift, image_translation, &
1710!$OMP ispin, jdir, target_atom) &
1711!$OMP SHARED(cell, composite_cross_force, composite_cross_image_virial, &
1712!$OMP composite_density_grad, composite_grad_grad, composite_grid_atom, &
1713!$OMP composite_grid_coords, composite_kin_grad, composite_nflat, cross_cutoff, &
1714!$OMP drho_h, drho_s, grid_atom, harmonics, image_shell, lsd, nspins, &
1715!$OMP particle_set, rho_h, rho_s, source_atom, tau_h, tau_s, &
1716!$OMP use_atom_composite_density, use_atom_composite_gradient, use_atom_composite_tau)
1717 ALLOCATE (composite_cross_force_local(3, SIZE(particle_set)))
1718 composite_cross_force_local = 0.0_dp
1719 composite_cross_image_virial_local = 0.0_dp
1720!$OMP DO SCHEDULE(STATIC)
1721 DO composite_row = 1, composite_nflat
1722 target_atom = composite_grid_atom(composite_row)
1723 IF (lsd) THEN
1724 cross_density_adjoint = 0.0_dp
1725 cross_grad_adjoint = 0.0_dp
1726 cross_kin_adjoint = 0.0_dp
1727 IF (use_atom_composite_density) THEN
1728 cross_density_adjoint(1:2) = &
1729 composite_density_grad(composite_row, 1:2)
1730 END IF
1731 IF (use_atom_composite_gradient) THEN
1732 cross_grad_adjoint(:, 1:2) = &
1733 composite_grad_grad(composite_row, :, 1:2)
1734 END IF
1735 IF (use_atom_composite_tau) THEN
1736 cross_kin_adjoint(1:2) = &
1737 composite_kin_grad(composite_row, 1:2)
1738 END IF
1739 ELSE
1740 cross_density_adjoint = 0.0_dp
1741 cross_grad_adjoint = 0.0_dp
1742 cross_kin_adjoint = 0.0_dp
1743 IF (use_atom_composite_density) THEN
1744 cross_density_adjoint(1) = 0.5_dp* &
1745 sum(composite_density_grad(composite_row, :))
1746 END IF
1747 IF (use_atom_composite_gradient) THEN
1748 DO idir = 1, 3
1749 cross_grad_adjoint(idir, 1) = 0.5_dp* &
1750 sum(composite_grad_grad(composite_row, idir, :))
1751 END DO
1752 END IF
1753 IF (use_atom_composite_tau) THEN
1754 cross_kin_adjoint(1) = 0.5_dp* &
1755 sum(composite_kin_grad(composite_row, :))
1756 END IF
1757 END IF
1758
1759 fractional = 0.0_dp
1760 DO idir = 1, 3
1761 DO jdir = 1, 3
1762 fractional(idir) = fractional(idir) + cell%h_inv(idir, jdir)* &
1763 (composite_grid_coords(jdir, composite_row) - &
1764 particle_set(source_atom)%r(jdir))
1765 END DO
1766 END DO
1767 CALL atom_grid_image_bounds(cell, fractional, cross_cutoff, image_shell, &
1768 image_lower, image_upper)
1769 DO image_i3 = image_lower(3), image_upper(3)
1770 DO image_i2 = image_lower(2), image_upper(2)
1771 DO image_i1 = image_lower(1), image_upper(1)
1772 image_shift = [image_i1, image_i2, image_i3]
1773 IF (target_atom == source_atom .AND. &
1774 all(image_shift == 0)) cycle
1775 image_translation = matmul( &
1776 cell%hmat, real(image_shift, dp))
1777 cross_displacement = &
1778 composite_grid_coords(:, composite_row) - &
1779 particle_set(source_atom)%r - image_translation
1780 IF (sqrt(sum(cross_displacement**2)) > cross_cutoff) cycle
1781 CALL interpolate_gapw_atom_grid_fields( &
1782 grid_atom, harmonics, cross_displacement, cross_cutoff, &
1783 nspins, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s, &
1784 cross_density, cross_grad, cross_kin, cross_density_spatial, &
1785 cross_grad_spatial, cross_kin_spatial)
1786 cross_spatial_derivative = 0.0_dp
1787 DO ispin = 1, nspins
1788 DO idir = 1, 3
1789 cross_spatial_derivative(idir) = &
1790 cross_spatial_derivative(idir) + &
1791 cross_density_adjoint(ispin)* &
1792 cross_density_spatial(idir, ispin) + &
1793 cross_kin_adjoint(ispin)* &
1794 cross_kin_spatial(idir, ispin)
1795 DO jdir = 1, 3
1796 cross_spatial_derivative(idir) = &
1797 cross_spatial_derivative(idir) + &
1798 cross_grad_adjoint(jdir, ispin)* &
1799 cross_grad_spatial(jdir, idir, ispin)
1800 END DO
1801 END DO
1802 END DO
1803 composite_cross_force_local(:, target_atom) = &
1804 composite_cross_force_local(:, target_atom) + &
1805 cross_spatial_derivative
1806 composite_cross_force_local(:, source_atom) = &
1807 composite_cross_force_local(:, source_atom) - &
1808 cross_spatial_derivative
1809 DO idir = 1, 3
1810 DO jdir = 1, 3
1811 composite_cross_image_virial_local(idir, jdir) = &
1812 composite_cross_image_virial_local(idir, jdir) + &
1813 cross_spatial_derivative(idir)*image_translation(jdir)
1814 END DO
1815 END DO
1816 END DO
1817 END DO
1818 END DO
1819 END DO
1820!$OMP END DO
1821!$OMP CRITICAL(skala_atom_composite_cross_reduction)
1822 composite_cross_force(:, :) = &
1823 composite_cross_force(:, :) + composite_cross_force_local(:, :)
1824 composite_cross_image_virial = composite_cross_image_virial + &
1825 composite_cross_image_virial_local
1826!$OMP END CRITICAL(skala_atom_composite_cross_reduction)
1827 DEALLOCATE (composite_cross_force_local)
1828!$OMP END PARALLEL
1829 END DO
1830 CALL release_tau_basis_cache(tau_basis_cache)
1831 END DO
1832 END IF
1833
1834 composite_explicit_force(:, :) = composite_model_atom_force(:, :) + &
1835 composite_grid_coord_force(:, :) + &
1836 composite_moving_smooth_force(:, :) + &
1837 composite_cross_force(:, :) + &
1838 composite_nlcc_center_force(:, :) + &
1839 composite_nlcc_target_force(:, :) + &
1840 composite_partition_force(:, :)
1841 ! CP2K stores +dE/dR in the electronic force components, while the
1842 ! virial is -dE/dstrain. Model coordinates, atom-grid centers,
1843 ! NLCC centers, and partition centers move affinely with their atoms.
1844 ! The smooth-field interpolation force is excluded here: its affine
1845 ! response is already in the PW stress and its non-affine local-grid
1846 ! correction is composite_interpolation_virial.
1847 DO iatom = 1, SIZE(particle_set)
1848 DO idir = 1, 3
1849 DO jdir = 1, 3
1850 composite_explicit_virial(idir, jdir) = &
1851 composite_explicit_virial(idir, jdir) - ( &
1852 composite_model_atom_force(idir, iatom) + &
1853 composite_grid_coord_force(idir, iatom) + &
1854 composite_cross_force(idir, iatom) + &
1855 composite_nlcc_center_force(idir, iatom) + &
1856 composite_nlcc_target_force(idir, iatom) + &
1857 composite_partition_force(idir, iatom))*particle_set(iatom)%r(jdir)
1858 END DO
1859 END DO
1860 END DO
1861 composite_explicit_virial = composite_explicit_virial + &
1862 composite_partition_strain_virial + &
1863 composite_cross_image_virial
1864 IF (ASSOCIATED(force)) THEN
1865 DO iatom = 1, SIZE(particle_set)
1866 ikind = composite_atom_kind(iatom)
1867 iat = composite_atom_kind_index(iatom)
1868 cpassert(ikind > 0 .AND. iat > 0)
1869 force(ikind)%rho_elec(:, iat) = force(ikind)%rho_elec(:, iat) + &
1870 composite_explicit_force(:, iatom)
1871 END DO
1872 END IF
1873 IF (use_virial) THEN
1874 ! The local radial vectors of the atom-centered quadrature do not
1875 ! deform with the periodic cell. Correct the standard PW response,
1876 ! which follows fixed fractional coordinates, by the corresponding
1877 ! non-affine interpolation derivative.
1878 virial%pv_xc = composite_feature_virial - composite_interpolation_virial
1879 virial%pv_gapw = virial%pv_gapw + composite_explicit_virial
1880 virial%pv_virial = virial%pv_virial + composite_explicit_virial
1881 END IF
1882 IF (native_grid_diagnostics) THEN
1883 CALL para_env%sum(composite_cross_force)
1884 CALL para_env%sum(composite_cross_image_virial)
1885 CALL para_env%sum(composite_model_atom_force)
1886 CALL para_env%sum(composite_grid_coord_force)
1887 CALL para_env%sum(composite_moving_smooth_force)
1888 CALL para_env%sum(composite_nlcc_center_force)
1889 CALL para_env%sum(composite_nlcc_target_force)
1890 CALL para_env%sum(composite_partition_force)
1891 CALL para_env%sum(composite_explicit_force)
1892 CALL para_env%sum(composite_explicit_virial)
1893 CALL para_env%sum(composite_feature_virial)
1894 CALL para_env%sum(composite_interpolation_virial)
1895 iw = cp_logger_get_default_io_unit()
1896 IF (iw > 0) THEN
1897 DO iatom = 1, SIZE(particle_set)
1898 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
1899 "SKALA_GPW| Atom-composite model-atom force", iatom, &
1900 composite_model_atom_force(:, iatom)
1901 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
1902 "SKALA_GPW| Atom-composite grid-coordinate force", iatom, &
1903 composite_grid_coord_force(:, iatom)
1904 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
1905 "SKALA_GPW| Atom-composite moving-smooth force", iatom, &
1906 composite_moving_smooth_force(:, iatom)
1907 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
1908 "SKALA_GPW| Atom-composite cross-region force", iatom, &
1909 composite_cross_force(:, iatom)
1910 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
1911 "SKALA_GPW| Atom-composite NLCC-center force", iatom, &
1912 composite_nlcc_center_force(:, iatom)
1913 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
1914 "SKALA_GPW| Atom-composite NLCC-target force", iatom, &
1915 composite_nlcc_target_force(:, iatom)
1916 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
1917 "SKALA_GPW| Atom-composite partition force", iatom, &
1918 composite_partition_force(:, iatom)
1919 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES20.12))") &
1920 "SKALA_GPW| Atom-composite explicit force", iatom, &
1921 composite_explicit_force(:, iatom)
1922 END DO
1923 WRITE (unit=iw, fmt="(T2,A)") &
1924 "SKALA_GPW| Atom-composite explicit virial"
1925 DO idir = 1, 3
1926 WRITE (unit=iw, fmt="(T2,A,1X,3ES20.10)") &
1927 "SKALA_GPW|", composite_explicit_virial(idir, :)
1928 END DO
1929 WRITE (unit=iw, fmt="(T2,A)") &
1930 "SKALA_GPW| Atom-composite cross-image virial"
1931 DO idir = 1, 3
1932 WRITE (unit=iw, fmt="(T2,A,1X,3ES20.10)") &
1933 "SKALA_GPW|", composite_cross_image_virial(idir, :)
1934 END DO
1935 WRITE (unit=iw, fmt="(T2,A)") &
1936 "SKALA_GPW| Atom-composite feature virial"
1937 DO idir = 1, 3
1938 WRITE (unit=iw, fmt="(T2,A,1X,3ES20.10)") &
1939 "SKALA_GPW|", composite_feature_virial(idir, :)
1940 END DO
1941 WRITE (unit=iw, fmt="(T2,A)") &
1942 "SKALA_GPW| Atom-composite interpolation virial"
1943 DO idir = 1, 3
1944 WRITE (unit=iw, fmt="(T2,A,1X,3ES20.10)") &
1945 "SKALA_GPW|", composite_interpolation_virial(idir, :)
1946 END DO
1947 END IF
1948 END IF
1949 DEALLOCATE (composite_atom_coord_grad, composite_atomic_grid_weight_grad, &
1950 composite_cross_force, &
1951 composite_explicit_force, composite_grid_coord_force, &
1952 composite_grid_coord_grad, composite_grid_weight_grad, &
1953 composite_model_atom_force, composite_moving_smooth_force, &
1954 composite_nlcc_center_force, &
1955 composite_nlcc_target_force, &
1956 composite_partition_datom, composite_partition_dstrain, &
1957 composite_partition_force, composite_partition_included)
1958 ELSE
1959 CALL skala_gapw_atom_composite_energy( &
1960 my_xc_section, para_env, composite_density, composite_grad, composite_kin, &
1961 composite_grid_coords, composite_grid_weights, composite_atomic_grid_weights, &
1962 composite_local_grid_sizes, composite_local_atom_coords, atom_composite_exc, &
1963 composite_density_grad, composite_grad_grad, composite_kin_grad)
1964 END IF
1965 composite_pw_nflat = product(smooth_rho_r(1)%pw_grid%npts)
1966 adjoint_nchannels = merge(2, 1, lsd)
1967 ALLOCATE (smooth_density_adjoint_storage(composite_pw_nflat, adjoint_nchannels), &
1968 smooth_grad_adjoint_storage(composite_pw_nflat, 3, adjoint_nchannels), &
1969 smooth_kin_adjoint_storage(composite_pw_nflat, adjoint_nchannels))
1970 smooth_density_adjoint_storage = 0.0_dp
1971 smooth_grad_adjoint_storage = 0.0_dp
1972 smooth_kin_adjoint_storage = 0.0_dp
1973 CALL build_native_grid_adjoint_bins( &
1974 smooth_rho_r(1)%pw_grid, cell, composite_grid_coords, &
1975 image_partition_atom_composite, adjoint_tile_count, adjoint_bin_offsets, &
1976 adjoint_bin_rows)
1977 adjoint_nbins = SIZE(adjoint_bin_offsets) - 1
1978!$OMP PARALLEL DO SCHEDULE(DYNAMIC) DEFAULT(NONE) &
1979!$OMP PRIVATE(adjoint_entry, composite_row, composite_smooth_density_adjoint_value, &
1980!$OMP composite_smooth_gradient_adjoint_value, composite_smooth_kin_adjoint_value, &
1981!$OMP adjoint_tile_lower, adjoint_tile_upper, interpolation_stencil) &
1982!$OMP SHARED(adjoint_bin_offsets, adjoint_bin_rows, adjoint_nbins, adjoint_nchannels, &
1983!$OMP adjoint_tile_count, cell, composite_density_grad, composite_grad_grad, &
1984!$OMP composite_grid_coords, composite_kin_grad, image_partition_atom_composite, lsd, &
1985!$OMP smooth_density_adjoint_storage, smooth_grad_adjoint_storage, &
1986!$OMP smooth_kin_adjoint_storage, smooth_rho_r, qs_env)
1987 DO adjoint_bin = 1, adjoint_nbins
1988 CALL native_grid_adjoint_tile_bounds( &
1989 smooth_rho_r(1)%pw_grid, adjoint_bin, adjoint_tile_count, &
1990 adjoint_tile_lower, adjoint_tile_upper)
1991 DO adjoint_entry = adjoint_bin_offsets(adjoint_bin), &
1992 adjoint_bin_offsets(adjoint_bin + 1) - 1
1993 composite_row = adjoint_bin_rows(adjoint_entry)
1994 IF (.NOT. fetch_native_grid_stencil(qs_env%native_grid_cache, composite_row, &
1995 composite_grid_coords(:, composite_row), interpolation_stencil)) THEN
1996 CALL create_native_grid_interpolation_stencil( &
1997 interpolation_stencil, smooth_rho_r(1)%pw_grid, cell, &
1998 composite_grid_coords(:, composite_row), image_partition_atom_composite)
1999 END IF
2000 IF (lsd) THEN
2001 composite_smooth_density_adjoint_value = &
2002 composite_density_grad(composite_row, :)
2003 composite_smooth_gradient_adjoint_value = &
2004 composite_grad_grad(composite_row, :, :)
2005 composite_smooth_kin_adjoint_value = composite_kin_grad(composite_row, :)
2006 ELSE
2007 composite_smooth_density_adjoint_value = 0.0_dp
2008 composite_smooth_gradient_adjoint_value = 0.0_dp
2009 composite_smooth_kin_adjoint_value = 0.0_dp
2010 composite_smooth_density_adjoint_value(1) = &
2011 sum(composite_density_grad(composite_row, :))
2012 composite_smooth_gradient_adjoint_value(:, 1) = &
2013 sum(composite_grad_grad(composite_row, :, :), dim=2)
2014 composite_smooth_kin_adjoint_value(1) = &
2015 sum(composite_kin_grad(composite_row, :))
2016 END IF
2017 CALL add_native_grid_fields_adjoint_tile( &
2018 smooth_density_adjoint_storage, smooth_grad_adjoint_storage, &
2019 smooth_kin_adjoint_storage, smooth_rho_r(1)%pw_grid, interpolation_stencil, &
2020 composite_smooth_density_adjoint_value, &
2021 composite_smooth_gradient_adjoint_value, &
2022 composite_smooth_kin_adjoint_value, adjoint_nchannels, &
2023 adjoint_tile_lower, adjoint_tile_upper)
2024 END DO
2025 END DO
2026!$OMP END PARALLEL DO
2027 DEALLOCATE (adjoint_bin_offsets, adjoint_bin_rows)
2028 smooth_input_contraction = 0.0_dp
2029!$OMP PARALLEL DO SCHEDULE(STATIC) REDUCTION(+:smooth_input_contraction) DEFAULT(NONE) &
2030!$OMP PRIVATE(composite_smooth_density_value, composite_smooth_gradient_value, &
2031!$OMP composite_smooth_kin_value, idir, ispin) &
2032!$OMP SHARED(composite_density_grad, composite_grad_grad, composite_kin_grad, composite_nflat, &
2033!$OMP composite_smooth_density_cache, composite_smooth_gradient_cache, &
2034!$OMP composite_smooth_kin_cache, lsd)
2035 DO composite_row = 1, composite_nflat
2036 composite_smooth_density_value = composite_smooth_density_cache(composite_row, :)
2037 composite_smooth_gradient_value = &
2038 composite_smooth_gradient_cache(composite_row, :, :)
2039 composite_smooth_kin_value = composite_smooth_kin_cache(composite_row, :)
2040 DO ispin = 1, 2
2041 IF (lsd) THEN
2042 smooth_input_contraction = smooth_input_contraction + &
2043 composite_density_grad(composite_row, ispin)* &
2044 composite_smooth_density_value(ispin) + &
2045 composite_kin_grad(composite_row, ispin)* &
2046 composite_smooth_kin_value(ispin)
2047 DO idir = 1, 3
2048 smooth_input_contraction = smooth_input_contraction + &
2049 composite_grad_grad(composite_row, idir, ispin)* &
2050 composite_smooth_gradient_value(idir, ispin)
2051 END DO
2052 ELSE
2053 smooth_input_contraction = smooth_input_contraction + 0.5_dp*( &
2054 composite_density_grad(composite_row, ispin)* &
2055 composite_smooth_density_value(1) + &
2056 composite_kin_grad(composite_row, ispin)* &
2057 composite_smooth_kin_value(1))
2058 DO idir = 1, 3
2059 smooth_input_contraction = smooth_input_contraction + 0.5_dp* &
2060 composite_grad_grad(composite_row, idir, ispin)* &
2061 composite_smooth_gradient_value(idir, 1)
2062 END DO
2063 END IF
2064 END DO
2065 END DO
2066!$OMP END PARALLEL DO
2067 smooth_density_adjoint => smooth_density_adjoint_storage
2068 smooth_grad_adjoint => smooth_grad_adjoint_storage
2069 smooth_kin_adjoint => smooth_kin_adjoint_storage
2070 ! The interpolation transpose still uses CP2K's global FFT-grid layout. Reduce only
2071 ! this PW adjoint; atom-grid feature rows and their model derivatives stay rank-local.
2072 CALL para_env%sum(smooth_density_adjoint)
2073 CALL para_env%sum(smooth_grad_adjoint)
2074 CALL para_env%sum(smooth_kin_adjoint)
2075 CALL para_env%sum(smooth_input_contraction)
2076 CALL build_vxc_from_feature_grads( &
2077 smooth_vxc_rho, smooth_vxc_tau, smooth_rho_r, auxbas_pw_pool, &
2078 smooth_density_adjoint, smooth_grad_adjoint, smooth_kin_adjoint, &
2079 xc_deriv_method_id, global_grid_layout=.true.)
2080 NULLIFY (smooth_density_adjoint, smooth_grad_adjoint, smooth_kin_adjoint)
2081 DEALLOCATE (smooth_density_adjoint_storage, smooth_grad_adjoint_storage, &
2082 smooth_kin_adjoint_storage)
2083 smooth_grid_contraction = 0.0_dp
2084 DO ispin = 1, nspins
2085 smooth_grid_contraction = smooth_grid_contraction + smooth_rho_r(1)%pw_grid%dvol* &
2086 (sum(smooth_vxc_rho(ispin)%array* &
2087 smooth_rho_r(ispin)%array) + &
2088 sum(smooth_vxc_tau(ispin)%array* &
2089 smooth_tau_r(ispin)%array))
2090 IF (atom_composite_reference) THEN
2091 cpassert(PRESENT(composite_vxc_rho))
2092 cpassert(PRESENT(composite_vxc_tau))
2093 cpassert(ASSOCIATED(composite_vxc_rho))
2094 cpassert(ASSOCIATED(composite_vxc_tau))
2095 cpassert(SIZE(composite_vxc_rho) == nspins)
2096 cpassert(SIZE(composite_vxc_tau) == nspins)
2097 CALL pw_axpy(smooth_vxc_rho(ispin), composite_vxc_rho(ispin), 1.0_dp)
2098 CALL pw_axpy(smooth_vxc_tau(ispin), composite_vxc_tau(ispin), 1.0_dp)
2099 END IF
2100 CALL auxbas_pw_pool%give_back_pw(smooth_vxc_rho(ispin))
2101 CALL auxbas_pw_pool%give_back_pw(smooth_vxc_tau(ispin))
2102 END DO
2103 CALL para_env%sum(smooth_grid_contraction)
2104 DEALLOCATE (smooth_vxc_rho, smooth_vxc_tau)
2105
2106 one_center_field_contraction = 0.0_dp
2107 one_center_matrix_contraction = 0.0_dp
2108 one_center_density_field_contraction = 0.0_dp
2109 one_center_density_matrix_contraction = 0.0_dp
2110 one_center_gradient_field_contraction = 0.0_dp
2111 one_center_gradient_matrix_contraction = 0.0_dp
2112 one_center_rho_grad_field_contraction = 0.0_dp
2113 one_center_rho_grad_matrix_contraction = 0.0_dp
2114 one_center_tau_field_contraction = 0.0_dp
2115 one_center_tau_matrix_contraction = 0.0_dp
2116 IF (.NOT. direct_valence_atom_composite) THEN
2117 DO ikind = 1, SIZE(atomic_kind_set)
2118 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
2119 NULLIFY (gth_potential, sgp_potential)
2120 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
2121 gth_potential=gth_potential, harmonics=harmonics, &
2122 grid_atom=grid_atom, sgp_potential=sgp_potential, &
2123 zatom=zatom, zeff=zeff)
2124 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
2125 IF (.NOT. native_skala_uses_one_center_kind( &
2126 paw_atom, gapw_representation, &
2127 ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential), &
2128 zeff, zatom)) cycle
2129
2130 nr = grid_atom%nr
2131 na = grid_atom%ng_sphere
2132 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
2133 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
2134 CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
2135 CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
2136 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
2137 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
2138 CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
2139 CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
2140 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
2141 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
2142 CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
2143 CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
2144 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
2145
2146 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
2147 DO iat = 1, natom
2148 iatom = atom_list(iat)
2149 source_matrix_local = iat >= bo(1) .AND. iat <= bo(2)
2150 rho_atom => my_rho_atom_set(iatom)
2151 NULLIFY (cpc_h, cpc_s, r_h, r_s, dr_h, dr_s, r_h_d, r_s_d, int_hh, int_ss)
2152 CALL get_rho_atom(rho_atom=rho_atom, cpc_h=cpc_h, cpc_s=cpc_s, &
2153 rho_rad_h=r_h, rho_rad_s=r_s, drho_rad_h=dr_h, &
2154 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, rho_rad_s_d=r_s_d, &
2155 ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
2156 rho_h = 0.0_dp
2157 rho_s = 0.0_dp
2158 drho_h = 0.0_dp
2159 drho_s = 0.0_dp
2160 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
2161 DO ir = 1, nr
2162 CALL calc_rho_angular(grid_atom, harmonics, nspins, .true., &
2163 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
2164 r_h_d, r_s_d, drho_h, drho_s)
2165 END DO
2166
2167 vxc_h = 0.0_dp
2168 vxc_s = 0.0_dp
2169 vxg_h = 0.0_dp
2170 vxg_s = 0.0_dp
2171 vtau_h = 0.0_dp
2172 vtau_s = 0.0_dp
2173 IF (source_matrix_local) THEN
2174 composite_row = composite_atom_start(iatom) - 1
2175 DO ir = 1, nr
2176 DO ia = 1, na
2177 composite_row = composite_row + 1
2178 IF (lsd) THEN
2179 IF (use_atom_composite_density) THEN
2180 vxc_h(ia, ir, 1:2) = composite_density_grad(composite_row, 1:2)
2181 vxc_s(ia, ir, 1:2) = composite_density_grad(composite_row, 1:2)
2182 END IF
2183 IF (use_atom_composite_gradient) THEN
2184 DO idir = 1, 3
2185 vxg_h(idir, ia, ir, 1:2) = &
2186 composite_grad_grad(composite_row, idir, 1:2)
2187 vxg_s(idir, ia, ir, 1:2) = &
2188 composite_grad_grad(composite_row, idir, 1:2)
2189 END DO
2190 END IF
2191 IF (use_atom_composite_tau) THEN
2192 vtau_h(ia, ir, 1:2) = composite_kin_grad(composite_row, 1:2)
2193 vtau_s(ia, ir, 1:2) = composite_kin_grad(composite_row, 1:2)
2194 END IF
2195 ELSE
2196 IF (use_atom_composite_density) THEN
2197 vxc_h(ia, ir, 1) = 0.5_dp* &
2198 sum(composite_density_grad(composite_row, :))
2199 vxc_s(ia, ir, 1) = vxc_h(ia, ir, 1)
2200 END IF
2201 IF (use_atom_composite_gradient) THEN
2202 DO idir = 1, 3
2203 vxg_h(idir, ia, ir, 1) = 0.5_dp* &
2204 sum(composite_grad_grad( &
2205 composite_row, idir, :))
2206 vxg_s(idir, ia, ir, 1) = vxg_h(idir, ia, ir, 1)
2207 END DO
2208 END IF
2209 IF (use_atom_composite_tau) THEN
2210 vtau_h(ia, ir, 1) = 0.5_dp* &
2211 sum(composite_kin_grad(composite_row, :))
2212 vtau_s(ia, ir, 1) = vtau_h(ia, ir, 1)
2213 END IF
2214 END IF
2215 END DO
2216 END DO
2217 cpassert(composite_row == composite_atom_end(iatom))
2218 END IF
2219
2220 cross_cutoff = gapw_atom_grid_support_radius( &
2221 grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s)
2222 IF (cross_cutoff > 0.0_dp) THEN
2223 image_shell = 0
2224 DO idir = 1, 3
2225 IF (cell%perd(idir) == 1) THEN
2226 image_shell(idir) = ceiling( &
2227 cross_cutoff*sqrt(sum(cell%h_inv(idir, :)**2))) + 1
2228 END IF
2229 END DO
2230 ! Each thread accumulates its own one-center potentials; combine them
2231 ! before the existing MPI reduction over target-atom owners.
2232!$OMP PARALLEL DEFAULT(NONE) &
2233!$OMP PRIVATE(image_lower, image_upper, composite_row, cross_density_adjoint, cross_displacement, &
2234!$OMP cross_grad_adjoint, cross_kin_adjoint, fractional, idir, image_i1, &
2235!$OMP image_i2, image_i3, image_shift, image_translation, jdir, target_atom, &
2236!$OMP vxc_h_local, vxc_s_local, vxg_h_local, vxg_s_local, vtau_h_local, vtau_s_local) &
2237!$OMP SHARED(cell, composite_density_grad, composite_grad_grad, composite_grid_atom, &
2238!$OMP composite_grid_coords, composite_kin_grad, composite_nflat, cross_cutoff, &
2239!$OMP grid_atom, harmonics, iatom, image_shell, lsd, nspins, particle_set, &
2240!$OMP use_atom_composite_density, use_atom_composite_gradient, use_atom_composite_tau, &
2241!$OMP vxc_h, vxc_s, vxg_h, vxg_s, vtau_h, vtau_s)
2242 ALLOCATE (vxc_h_local, mold=vxc_h)
2243 ALLOCATE (vxc_s_local, mold=vxc_s)
2244 ALLOCATE (vxg_h_local, mold=vxg_h)
2245 ALLOCATE (vxg_s_local, mold=vxg_s)
2246 ALLOCATE (vtau_h_local, mold=vtau_h)
2247 ALLOCATE (vtau_s_local, mold=vtau_s)
2248 vxc_h_local = 0.0_dp
2249 vxc_s_local = 0.0_dp
2250 vxg_h_local = 0.0_dp
2251 vxg_s_local = 0.0_dp
2252 vtau_h_local = 0.0_dp
2253 vtau_s_local = 0.0_dp
2254!$OMP DO SCHEDULE(STATIC)
2255 DO composite_row = 1, composite_nflat
2256 target_atom = composite_grid_atom(composite_row)
2257 cross_density_adjoint = 0.0_dp
2258 cross_grad_adjoint = 0.0_dp
2259 cross_kin_adjoint = 0.0_dp
2260 IF (lsd) THEN
2261 IF (use_atom_composite_density) THEN
2262 cross_density_adjoint(1:2) = &
2263 composite_density_grad(composite_row, 1:2)
2264 END IF
2265 IF (use_atom_composite_gradient) THEN
2266 cross_grad_adjoint(:, 1:2) = &
2267 composite_grad_grad(composite_row, :, 1:2)
2268 END IF
2269 IF (use_atom_composite_tau) THEN
2270 cross_kin_adjoint(1:2) = &
2271 composite_kin_grad(composite_row, 1:2)
2272 END IF
2273 ELSE
2274 IF (use_atom_composite_density) THEN
2275 cross_density_adjoint(1) = 0.5_dp* &
2276 sum(composite_density_grad(composite_row, :))
2277 END IF
2278 IF (use_atom_composite_gradient) THEN
2279 DO idir = 1, 3
2280 cross_grad_adjoint(idir, 1) = 0.5_dp* &
2281 sum(composite_grad_grad(composite_row, idir, :))
2282 END DO
2283 END IF
2284 IF (use_atom_composite_tau) THEN
2285 cross_kin_adjoint(1) = 0.5_dp* &
2286 sum(composite_kin_grad(composite_row, :))
2287 END IF
2288 END IF
2289 fractional = 0.0_dp
2290 DO idir = 1, 3
2291 DO jdir = 1, 3
2292 fractional(idir) = fractional(idir) + cell%h_inv(idir, jdir)* &
2293 (composite_grid_coords(jdir, composite_row) - &
2294 particle_set(iatom)%r(jdir))
2295 END DO
2296 END DO
2297 CALL atom_grid_image_bounds(cell, fractional, cross_cutoff, image_shell, &
2298 image_lower, image_upper)
2299 DO image_i3 = image_lower(3), image_upper(3)
2300 DO image_i2 = image_lower(2), image_upper(2)
2301 DO image_i1 = image_lower(1), image_upper(1)
2302 image_shift = [image_i1, image_i2, image_i3]
2303 IF (target_atom == iatom .AND. all(image_shift == 0)) cycle
2304 image_translation = matmul( &
2305 cell%hmat, real(image_shift, dp))
2306 cross_displacement = &
2307 composite_grid_coords(:, composite_row) - &
2308 particle_set(iatom)%r - image_translation
2309 ! Reject inactive images before allocating interpolation scratch.
2310 IF (sqrt(sum(cross_displacement**2)) > cross_cutoff) cycle
2311 CALL add_gapw_atom_grid_interpolation_adjoint( &
2312 grid_atom, harmonics, cross_displacement, cross_cutoff, &
2313 nspins, cross_density_adjoint, cross_grad_adjoint, &
2314 cross_kin_adjoint, vxc_h_local, vxc_s_local, vxg_h_local, vxg_s_local, &
2315 vtau_h_local, vtau_s_local)
2316 END DO
2317 END DO
2318 END DO
2319 END DO
2320!$OMP END DO
2321!$OMP CRITICAL(skala_atom_composite_adjoint_reduction)
2322 vxc_h = vxc_h + vxc_h_local
2323 vxc_s = vxc_s + vxc_s_local
2324 vxg_h = vxg_h + vxg_h_local
2325 vxg_s = vxg_s + vxg_s_local
2326 vtau_h = vtau_h + vtau_h_local
2327 vtau_s = vtau_s + vtau_s_local
2328!$OMP END CRITICAL(skala_atom_composite_adjoint_reduction)
2329 DEALLOCATE (vxc_h_local, vxc_s_local, vxg_h_local, vxg_s_local, vtau_h_local, vtau_s_local)
2330!$OMP END PARALLEL
2331 END IF
2332
2333 ! Model rows are distributed by target atom, while CP2K stores each
2334 ! one-center matrix on the rank owning its source atom. Sum the exact
2335 ! interpolation transpose before forming that matrix.
2336 CALL para_env%sum(vxc_h)
2337 CALL para_env%sum(vxc_s)
2338 CALL para_env%sum(vxg_h)
2339 CALL para_env%sum(vxg_s)
2340 CALL para_env%sum(vtau_h)
2341 CALL para_env%sum(vtau_s)
2342 IF (.NOT. source_matrix_local) cycle
2343
2344 one_center_rho_grad_field_contraction = &
2345 one_center_rho_grad_field_contraction + sum(vxc_h*(rho_h - rho_s)) + &
2346 sum(vxg_h*(drho_h(1:3, :, :, :) - drho_s(1:3, :, :, :)))
2347 one_center_density_field_contraction = one_center_density_field_contraction + &
2348 sum(vxc_h*(rho_h - rho_s))
2349 one_center_gradient_field_contraction = one_center_gradient_field_contraction + &
2350 sum(vxg_h*(drho_h(1:3, :, :, :) - &
2351 drho_s(1:3, :, :, :)))
2352 one_center_tau_field_contraction = one_center_tau_field_contraction + &
2353 sum(vtau_h*(tau_h - tau_s))
2354
2355 ALLOCATE (composite_int_h(SIZE(int_hh(1)%r_coef, 1), &
2356 SIZE(int_hh(1)%r_coef, 2), nspins), &
2357 composite_int_s(SIZE(int_ss(1)%r_coef, 1), &
2358 SIZE(int_ss(1)%r_coef, 2), nspins))
2359 DO ispin = 1, nspins
2360 composite_int_h(:, :, ispin) = int_hh(ispin)%r_coef
2361 composite_int_s(:, :, ispin) = int_ss(ispin)%r_coef
2362 int_hh(ispin)%r_coef = 0.0_dp
2363 int_ss(ispin)%r_coef = 0.0_dp
2364 END DO
2365 CALL gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, &
2366 grid_atom, basis_1c, harmonics, nspins)
2367 DO ispin = 1, nspins
2368 one_center_density_matrix_contraction = &
2369 one_center_density_matrix_contraction + &
2370 contract_one_center_matrix(cpc_h(ispin)%r_coef, &
2371 int_hh(ispin)%r_coef, &
2372 tau_basis_cache%n2oindex) - &
2373 contract_one_center_matrix(cpc_s(ispin)%r_coef, &
2374 int_ss(ispin)%r_coef, &
2375 tau_basis_cache%n2oindex)
2376 int_hh(ispin)%r_coef = 0.0_dp
2377 int_ss(ispin)%r_coef = 0.0_dp
2378 END DO
2379 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
2380 grid_atom, basis_1c, harmonics, nspins)
2381 DO ispin = 1, nspins
2382 one_center_rho_grad_matrix_contraction = &
2383 one_center_rho_grad_matrix_contraction + &
2384 contract_one_center_matrix(cpc_h(ispin)%r_coef, &
2385 int_hh(ispin)%r_coef, &
2386 tau_basis_cache%n2oindex) - &
2387 contract_one_center_matrix(cpc_s(ispin)%r_coef, &
2388 int_ss(ispin)%r_coef, &
2389 tau_basis_cache%n2oindex)
2390 END DO
2391 CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, tau_basis_cache, nspins)
2392 DO ispin = 1, nspins
2393 one_center_matrix_contraction = one_center_matrix_contraction + &
2394 contract_one_center_matrix(cpc_h(ispin)%r_coef, &
2395 int_hh(ispin)%r_coef, &
2396 tau_basis_cache%n2oindex) - &
2397 contract_one_center_matrix(cpc_s(ispin)%r_coef, &
2398 int_ss(ispin)%r_coef, &
2399 tau_basis_cache%n2oindex)
2400 IF (.NOT. atom_composite_reference) THEN
2401 int_hh(ispin)%r_coef = composite_int_h(:, :, ispin)
2402 int_ss(ispin)%r_coef = composite_int_s(:, :, ispin)
2403 END IF
2404 END DO
2405 DEALLOCATE (composite_int_h, composite_int_s)
2406 END DO
2407 CALL release_tau_basis_cache(tau_basis_cache)
2408 END DO
2409 END IF
2410 CALL para_env%sum(one_center_density_field_contraction)
2411 CALL para_env%sum(one_center_density_matrix_contraction)
2412 CALL para_env%sum(one_center_gradient_field_contraction)
2413 CALL para_env%sum(one_center_matrix_contraction)
2414 CALL para_env%sum(one_center_rho_grad_field_contraction)
2415 CALL para_env%sum(one_center_rho_grad_matrix_contraction)
2416 CALL para_env%sum(one_center_tau_field_contraction)
2417 one_center_field_contraction = one_center_rho_grad_field_contraction + &
2418 one_center_tau_field_contraction
2419 one_center_gradient_matrix_contraction = one_center_rho_grad_matrix_contraction - &
2420 one_center_density_matrix_contraction
2421 one_center_tau_matrix_contraction = one_center_matrix_contraction - &
2422 one_center_rho_grad_matrix_contraction
2423 feature_component_analytic(1) = sum(composite_density_grad*composite_density)
2424 feature_component_analytic(2) = sum(composite_grad_grad*composite_grad)
2425 feature_component_analytic(3) = sum(composite_kin_grad*composite_kin)
2426 feature_component_analytic(4) = sum(feature_component_analytic(1:3))
2427 CALL para_env%sum(feature_component_analytic)
2428 one_center_tensor_contraction = feature_component_analytic(4) - &
2429 smooth_input_contraction
2430 IF (atom_composite_diagnostic) THEN
2431 DO icomponent = 1, 4
2432 SELECT CASE (icomponent)
2433 CASE (1)
2434 composite_density = (1.0_dp + feature_vxc_step)*composite_density
2435 CASE (2)
2436 composite_grad = (1.0_dp + feature_vxc_step)*composite_grad
2437 CASE (3)
2438 composite_kin = (1.0_dp + feature_vxc_step)*composite_kin
2439 CASE (4)
2440 composite_density = (1.0_dp + feature_vxc_step)*composite_density
2441 composite_grad = (1.0_dp + feature_vxc_step)*composite_grad
2442 composite_kin = (1.0_dp + feature_vxc_step)*composite_kin
2443 END SELECT
2444 CALL skala_gapw_atom_composite_energy( &
2445 my_xc_section, para_env, composite_density, composite_grad, composite_kin, &
2446 composite_grid_coords, composite_grid_weights, composite_atomic_grid_weights, &
2447 composite_local_grid_sizes, composite_local_atom_coords, feature_vxc_plus)
2448 SELECT CASE (icomponent)
2449 CASE (1)
2450 composite_density = ((1.0_dp - feature_vxc_step)/ &
2451 (1.0_dp + feature_vxc_step))*composite_density
2452 CASE (2)
2453 composite_grad = ((1.0_dp - feature_vxc_step)/ &
2454 (1.0_dp + feature_vxc_step))*composite_grad
2455 CASE (3)
2456 composite_kin = ((1.0_dp - feature_vxc_step)/ &
2457 (1.0_dp + feature_vxc_step))*composite_kin
2458 CASE (4)
2459 composite_density = ((1.0_dp - feature_vxc_step)/ &
2460 (1.0_dp + feature_vxc_step))*composite_density
2461 composite_grad = ((1.0_dp - feature_vxc_step)/ &
2462 (1.0_dp + feature_vxc_step))*composite_grad
2463 composite_kin = ((1.0_dp - feature_vxc_step)/ &
2464 (1.0_dp + feature_vxc_step))*composite_kin
2465 END SELECT
2466 CALL skala_gapw_atom_composite_energy( &
2467 my_xc_section, para_env, composite_density, composite_grad, composite_kin, &
2468 composite_grid_coords, composite_grid_weights, composite_atomic_grid_weights, &
2469 composite_local_grid_sizes, composite_local_atom_coords, feature_vxc_minus)
2470 SELECT CASE (icomponent)
2471 CASE (1)
2472 composite_density = composite_density/(1.0_dp - feature_vxc_step)
2473 CASE (2)
2474 composite_grad = composite_grad/(1.0_dp - feature_vxc_step)
2475 CASE (3)
2476 composite_kin = composite_kin/(1.0_dp - feature_vxc_step)
2477 CASE (4)
2478 composite_density = composite_density/(1.0_dp - feature_vxc_step)
2479 composite_grad = composite_grad/(1.0_dp - feature_vxc_step)
2480 composite_kin = composite_kin/(1.0_dp - feature_vxc_step)
2481 END SELECT
2482 feature_component_fd(icomponent) = &
2483 (feature_vxc_plus - feature_vxc_minus)/(2.0_dp*feature_vxc_step)
2484 END DO
2485 feature_vxc_analytic = feature_component_analytic(4)
2486 feature_vxc_fd = feature_component_fd(4)
2487 iw = cp_logger_get_default_io_unit()
2488 IF (iw > 0) THEN
2489 WRITE (unit=iw, fmt="(T2,A,1X,I0)") &
2490 "SKALA_GPW| Atom-composite reference components", atom_composite_components
2491 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2492 "SKALA_GPW| Atom-composite reference electrons", atom_composite_nelec
2493 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2494 "SKALA_GPW| Atom-composite reference XC energy", atom_composite_exc
2495 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2496 "SKALA_GPW| Atom-composite feature VXC contraction", feature_vxc_analytic
2497 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2498 "SKALA_GPW| Atom-composite feature finite difference", feature_vxc_fd
2499 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2500 "SKALA_GPW| Atom-composite feature VXC error", &
2501 feature_vxc_analytic - feature_vxc_fd
2502 WRITE (unit=iw, fmt="(T2,A,3(1X,ES24.16))") &
2503 "SKALA_GPW| Atom-composite smooth PW adjoint", smooth_input_contraction, &
2504 smooth_grid_contraction, smooth_grid_contraction - smooth_input_contraction
2505 WRITE (unit=iw, fmt="(T2,A,4(1X,ES24.16))") &
2506 "SKALA_GPW| Atom-composite one-center adjoint", one_center_tensor_contraction, &
2507 one_center_field_contraction, one_center_matrix_contraction, &
2508 one_center_matrix_contraction - one_center_tensor_contraction
2509 WRITE (unit=iw, fmt="(T2,A,4(1X,ES24.16))") &
2510 "SKALA_GPW| Atom-composite one-center channels", &
2511 one_center_rho_grad_field_contraction, one_center_rho_grad_matrix_contraction, &
2512 one_center_tau_field_contraction, one_center_tau_matrix_contraction
2513 WRITE (unit=iw, fmt="(T2,A,4(1X,ES24.16))") &
2514 "SKALA_GPW| Atom-composite one-center rho-gradient", &
2515 one_center_density_field_contraction, one_center_density_matrix_contraction, &
2516 one_center_gradient_field_contraction, one_center_gradient_matrix_contraction
2517 DO icomponent = 1, 3
2518 WRITE (unit=iw, fmt="(T2,A,1X,I0,3(1X,ES24.16))") &
2519 "SKALA_GPW| Atom-composite component VXC", icomponent, &
2520 feature_component_analytic(icomponent), feature_component_fd(icomponent), &
2521 feature_component_analytic(icomponent) - feature_component_fd(icomponent)
2522 END DO
2523 END IF
2524 END IF
2525 IF (atom_composite_reference) THEN
2526 exc1 = atom_composite_exc
2527 IF (native_grid_diagnostics) THEN
2528 IF (composite_nflat > 0) THEN
2529 composite_density_min = minval(composite_density)
2530 composite_density_max = maxval(composite_density)
2531 composite_kin_min = minval(composite_kin)
2532 composite_kin_max = maxval(composite_kin)
2533 composite_grad_max = maxval(abs(composite_grad))
2534 ELSE
2535 composite_density_min = huge(1.0_dp)
2536 composite_density_max = -huge(1.0_dp)
2537 composite_kin_min = huge(1.0_dp)
2538 composite_kin_max = -huge(1.0_dp)
2539 composite_grad_max = 0.0_dp
2540 END IF
2541 composite_tau_integral = &
2542 sum(composite_grid_weights*sum(composite_kin, dim=2))
2543 CALL para_env%min(composite_density_min)
2544 CALL para_env%max(composite_density_max)
2545 CALL para_env%min(composite_kin_min)
2546 CALL para_env%max(composite_kin_max)
2547 CALL para_env%max(composite_grad_max)
2548 CALL para_env%sum(composite_tau_integral)
2549 END IF
2550 iw = cp_logger_get_default_io_unit()
2551 IF (iw > 0) THEN
2552 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2553 "SKALA_GPW| Active atom-composite XC energy", atom_composite_exc
2554 IF (native_grid_diagnostics) THEN
2555 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2556 "SKALA_GPW| Active atom-composite electrons", atom_composite_nelec
2557 WRITE (unit=iw, fmt="(T2,A,2(1X,ES24.16))") &
2558 "SKALA_GPW| Active atom-composite density range", &
2559 composite_density_min, composite_density_max
2560 WRITE (unit=iw, fmt="(T2,A,2(1X,ES24.16))") &
2561 "SKALA_GPW| Active atom-composite tau range", &
2562 composite_kin_min, composite_kin_max
2563 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2564 "SKALA_GPW| Active atom-composite tau integral", &
2565 composite_tau_integral
2566 WRITE (unit=iw, fmt="(T2,A,1X,ES24.16)") &
2567 "SKALA_GPW| Active atom-composite max gradient", &
2568 composite_grad_max
2569 END IF
2570 END IF
2571 END IF
2572 CALL xc_rho_set_release(smooth_rho_set, pw_pool=auxbas_pw_pool)
2573 DEALLOCATE (composite_atomic_grid_sizes, composite_atom_kind, &
2574 composite_atom_kind_index, composite_atom_start, composite_atom_end, &
2575 composite_atom_coords, composite_local_atoms, composite_local_grid_sizes, &
2576 composite_local_atom_coords, composite_grid_atom, &
2577 composite_partition_weights, composite_partition_atom_coords, &
2578 composite_distances, composite_density, composite_grad, composite_kin, &
2579 composite_smooth_density_cache, composite_smooth_gradient_cache, &
2580 composite_smooth_kin_cache, &
2581 composite_density_grad, composite_grad_grad, composite_kin_grad, &
2582 composite_grid_coords, composite_grid_weights, &
2583 composite_base_grid_weights, composite_atomic_grid_weights)
2584 IF (lsd) THEN
2585 DEALLOCATE (composite_smooth_rhoa, composite_smooth_rhob, &
2586 composite_smooth_tau_a, composite_smooth_tau_b)
2587 DO idir = 1, 3
2588 DEALLOCATE (composite_smooth_drhoa(idir)%array, &
2589 composite_smooth_drhob(idir)%array)
2590 END DO
2591 ELSE
2592 DEALLOCATE (composite_smooth_rho, composite_smooth_tau)
2593 DO idir = 1, 3
2594 DEALLOCATE (composite_smooth_drho(idir)%array)
2595 END DO
2596 END IF
2597 END IF
2598
2599 IF (.NOT. atom_composite_reference) CALL para_env%sum(exc1)
2600
2601 IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
2602 IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
2603 IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
2604 IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
2605
2606 IF (gradient_f) THEN
2607 IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
2608 IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
2609 IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
2610 IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
2611 END IF
2612
2613 IF (tau_f) THEN
2614 IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
2615 IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
2616 IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
2617 IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
2618 END IF
2619
2620 END IF !xc_none
2621
2622 CALL timestop(handle)
2623
2624 END SUBROUTINE calculate_vxc_atom
2625
2626! **************************************************************************************************
2627!> \brief Add the GAPW one-center correction to CDFT values and operators.
2628!> \param qs_env Quickstep environment
2629!> \param energy_only skip construction of the CDFT one-center operator
2630!> \param calculate_forces evaluate explicit derivatives of the partition weights
2631!> \param values constraint values from the hard-minus-soft one-center densities
2632!> \param electronic_charge optional one-center corrections to atomic populations
2633!> \param operator_group optional group for which to build the unscaled weight operator
2634!> \param rho_atom_operator_set optional destination for the one-center operator integrals
2635! **************************************************************************************************
2636 SUBROUTINE gapw_cdft_one_center(qs_env, energy_only, calculate_forces, values, &
2637 electronic_charge, operator_group, rho_atom_operator_set)
2638 TYPE(qs_environment_type), POINTER :: qs_env
2639 LOGICAL, INTENT(IN) :: energy_only, calculate_forces
2640 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: values
2641 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT), &
2642 OPTIONAL :: electronic_charge
2643 INTEGER, INTENT(IN), OPTIONAL :: operator_group
2644 TYPE(rho_atom_type), DIMENSION(:), POINTER, &
2645 OPTIONAL :: rho_atom_operator_set
2646
2647 INTEGER :: atom, channel, ia, iat, igroup, ikind, &
2648 ir, natom, natom_kind, nspins
2649 INTEGER, DIMENSION(2) :: atom_bounds
2650 INTEGER, DIMENSION(:), POINTER :: atom_list
2651 LOGICAL :: lsd, paw_atom
2652 REAL(kind=dp) :: delta_density, point_factor, spin_factor
2653 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: atomic_weights, group_weights
2654 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: explicit_derivative, &
2655 group_point_derivative
2656 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: group_atom_derivative
2657 REAL(kind=dp), DIMENSION(3) :: point
2658 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s
2659 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s
2660 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: vlocal_h, vlocal_s
2661 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
2662 TYPE(cdft_control_type), POINTER :: cdft_control
2663 TYPE(cdft_point_context_type) :: context
2664 TYPE(dft_control_type), POINTER :: dft_control
2665 TYPE(grid_atom_type), POINTER :: grid_atom
2666 TYPE(gto_basis_set_type), POINTER :: basis_1c
2667 TYPE(harmonics_atom_type), POINTER :: harmonics
2668 TYPE(mp_para_env_type), POINTER :: para_env
2669 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
2670 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
2671 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
2672 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
2673 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
2674 TYPE(rho_atom_type), DIMENSION(:), POINTER :: operator_atom_set, rho_atom_set
2675 TYPE(rho_atom_type), POINTER :: operator_atom, rho_atom
2676
2677 NULLIFY (atom_list, atomic_kind_set, basis_1c, cdft_control, dft_control, force, &
2678 grid_atom, harmonics, int_hh, int_ss, para_env, particle_set, r_h, r_s, &
2679 dr_h, dr_s, r_h_d, r_s_d, rho_h, rho_s, drho_h, drho_s, &
2680 operator_atom, operator_atom_set, rho_atom, rho_atom_set, qs_kind_set, &
2681 vlocal_h, vlocal_s)
2682 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, dft_control=dft_control, &
2683 force=force, natom=natom, para_env=para_env, particle_set=particle_set, &
2684 qs_kind_set=qs_kind_set, rho_atom_set=rho_atom_set)
2685 cpassert(ASSOCIATED(atomic_kind_set))
2686 cpassert(ASSOCIATED(dft_control))
2687 cpassert(ASSOCIATED(para_env))
2688 cpassert(ASSOCIATED(particle_set))
2689 cpassert(ASSOCIATED(qs_kind_set))
2690 cpassert(ASSOCIATED(rho_atom_set))
2691 operator_atom_set => rho_atom_set
2692 IF (PRESENT(rho_atom_operator_set)) operator_atom_set => rho_atom_operator_set
2693 cpassert(ASSOCIATED(operator_atom_set))
2694 cdft_control => dft_control%qs_control%cdft_control
2695 cpassert(ASSOCIATED(cdft_control))
2696 nspins = dft_control%nspins
2697 lsd = dft_control%lsd
2698 cpassert(SIZE(values) == SIZE(cdft_control%group))
2699 IF (PRESENT(operator_group)) THEN
2700 cpassert(operator_group >= 1 .AND. operator_group <= SIZE(cdft_control%group))
2701 END IF
2702 IF (PRESENT(electronic_charge)) THEN
2703 cpassert(SIZE(electronic_charge, 1) == natom)
2704 cpassert(SIZE(electronic_charge, 2) == nspins)
2705 END IF
2706
2707 CALL cdft_point_context_create(qs_env, context, calculate_forces)
2708 ALLOCATE (group_weights(context%ngroup), group_point_derivative(3, context%ngroup), &
2709 group_atom_derivative(3, natom, context%ngroup), &
2710 explicit_derivative(3, natom))
2711 IF (PRESENT(electronic_charge)) ALLOCATE (atomic_weights(natom))
2712 values = 0.0_dp
2713 explicit_derivative = 0.0_dp
2714 IF (PRESENT(electronic_charge)) electronic_charge = 0.0_dp
2715
2716 DO ikind = 1, SIZE(atomic_kind_set)
2717 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom_kind)
2718 CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, grid_atom=grid_atom, &
2719 harmonics=harmonics)
2720 IF (.NOT. paw_atom) cycle
2721 CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
2722 cpassert(ASSOCIATED(grid_atom))
2723 cpassert(ASSOCIATED(harmonics))
2724 cpassert(ASSOCIATED(basis_1c))
2725 ALLOCATE (vlocal_h(grid_atom%ng_sphere, grid_atom%nr, nspins), &
2726 vlocal_s(grid_atom%ng_sphere, grid_atom%nr, nspins))
2727
2728 atom_bounds = get_limit(natom_kind, para_env%num_pe, para_env%mepos)
2729 DO iat = atom_bounds(1), atom_bounds(2)
2730 atom = atom_list(iat)
2731 rho_atom => rho_atom_set(atom)
2732 NULLIFY (r_h, r_s)
2733 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
2734 CALL reallocate(rho_h, 1, grid_atom%ng_sphere, 1, grid_atom%nr, 1, nspins)
2735 CALL reallocate(rho_s, 1, grid_atom%ng_sphere, 1, grid_atom%nr, 1, nspins)
2736 rho_h = 0.0_dp
2737 rho_s = 0.0_dp
2738 DO ir = 1, grid_atom%nr
2739 CALL calc_rho_angular(grid_atom, harmonics, nspins, .false., ir, r_h, r_s, &
2740 rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
2741 END DO
2742 vlocal_h = 0.0_dp
2743 vlocal_s = 0.0_dp
2744
2745 DO ir = 1, grid_atom%nr
2746 DO ia = 1, grid_atom%ng_sphere
2747 point(1) = particle_set(atom)%r(1) + grid_atom%rad(ir)* &
2748 grid_atom%sin_pol(ia)*grid_atom%cos_azi(ia)
2749 point(2) = particle_set(atom)%r(2) + grid_atom%rad(ir)* &
2750 grid_atom%sin_pol(ia)*grid_atom%sin_azi(ia)
2751 point(3) = particle_set(atom)%r(3) + grid_atom%rad(ir)*grid_atom%cos_pol(ia)
2752 IF (PRESENT(electronic_charge)) THEN
2753 CALL cdft_point_weights(context, point, group_weights, &
2754 group_point_derivative, group_atom_derivative, &
2755 atomic_weights)
2756 ELSE
2757 CALL cdft_point_weights(context, point, group_weights, &
2758 group_point_derivative, group_atom_derivative)
2759 END IF
2760 DO channel = 1, nspins
2761 delta_density = rho_h(ia, ir, channel) - rho_s(ia, ir, channel)
2762 IF (PRESENT(electronic_charge)) THEN
2763 electronic_charge(:, channel) = electronic_charge(:, channel) + &
2764 grid_atom%weight(ia, ir)*atomic_weights* &
2765 delta_density
2766 END IF
2767 point_factor = 0.0_dp
2768 DO igroup = 1, context%ngroup
2769 spin_factor = cdft_spin_factor( &
2770 cdft_control%group(igroup)%constraint_type, channel, lsd)
2771 IF (PRESENT(operator_group)) THEN
2772 IF (igroup == operator_group) THEN
2773 point_factor = point_factor + group_weights(igroup)
2774 END IF
2775 ELSE
2776 point_factor = point_factor + cdft_control%strength(igroup)* &
2777 group_weights(igroup)*spin_factor
2778 END IF
2779 values(igroup) = values(igroup) + grid_atom%weight(ia, ir)* &
2780 group_weights(igroup)*delta_density*spin_factor
2781 IF (calculate_forces) THEN
2782 explicit_derivative(:, :) = &
2783 explicit_derivative + grid_atom%weight(ia, ir)* &
2784 cdft_control%strength(igroup)*delta_density* &
2785 spin_factor*group_atom_derivative(:, :, igroup)
2786 explicit_derivative(:, atom) = explicit_derivative(:, atom) + &
2787 grid_atom%weight(ia, ir)* &
2788 cdft_control%strength(igroup)*delta_density* &
2789 spin_factor*group_point_derivative(:, igroup)
2790 END IF
2791 END DO
2792 vlocal_h(ia, ir, channel) = grid_atom%weight(ia, ir)*point_factor
2793 vlocal_s(ia, ir, channel) = vlocal_h(ia, ir, channel)
2794 END DO
2795 END DO
2796 END DO
2797
2798 IF (.NOT. energy_only) THEN
2799 NULLIFY (int_hh, int_ss)
2800 operator_atom => operator_atom_set(atom)
2801 CALL get_rho_atom(rho_atom=operator_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
2802 CALL gavxcgb_nogc(vlocal_h, vlocal_s, int_hh, int_ss, grid_atom, &
2803 basis_1c, harmonics, nspins)
2804 END IF
2805 END DO
2806 DEALLOCATE (vlocal_h, vlocal_s)
2807 END DO
2808
2809 CALL para_env%sum(values)
2810 CALL para_env%sum(explicit_derivative)
2811 IF (PRESENT(electronic_charge)) CALL para_env%sum(electronic_charge)
2812 IF (calculate_forces .AND. ASSOCIATED(force) .AND. para_env%is_source()) THEN
2813 DO ikind = 1, SIZE(atomic_kind_set)
2814 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom_kind)
2815 DO iat = 1, natom_kind
2816 atom = atom_list(iat)
2817 force(ikind)%rho_elec(:, iat) = force(ikind)%rho_elec(:, iat) + &
2818 explicit_derivative(:, atom)
2819 END DO
2820 END DO
2821 END IF
2822
2823 IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
2824 IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
2825 IF (ALLOCATED(atomic_weights)) DEALLOCATE (atomic_weights)
2826 DEALLOCATE (explicit_derivative, group_atom_derivative, group_point_derivative, group_weights)
2827 CALL cdft_point_context_release(context)
2828
2829 CONTAINS
2830
2831! **************************************************************************************************
2832!> \brief ...
2833!> \param constraint_type ...
2834!> \param channel ...
2835!> \param lsd ...
2836!> \return ...
2837! **************************************************************************************************
2838 FUNCTION cdft_spin_factor(constraint_type, channel, lsd) RESULT(factor)
2839 INTEGER, INTENT(IN) :: constraint_type, channel
2840 LOGICAL, INTENT(IN) :: lsd
2841 REAL(kind=dp) :: factor
2842
2843 SELECT CASE (constraint_type)
2844 CASE (cdft_charge_constraint)
2845 factor = 1.0_dp
2846 CASE (cdft_magnetization_constraint)
2847 cpassert(lsd)
2848 factor = merge(1.0_dp, -1.0_dp, channel == 1)
2849 CASE (cdft_alpha_constraint)
2850 cpassert(lsd)
2851 factor = merge(1.0_dp, 0.0_dp, channel == 1)
2852 CASE (cdft_beta_constraint)
2853 cpassert(lsd)
2854 factor = merge(1.0_dp, 0.0_dp, channel == 2)
2855 CASE DEFAULT
2856 cpabort("Unknown CDFT constraint type.")
2857 END SELECT
2858 END FUNCTION cdft_spin_factor
2859
2860 END SUBROUTINE gapw_cdft_one_center
2861
2862! **************************************************************************************************
2863!> \brief Contract a compact one-center density matrix with an integral in the padded old basis.
2864!> \param density_matrix ...
2865!> \param integral_matrix ...
2866!> \param new_to_old ...
2867!> \return ...
2868! **************************************************************************************************
2869 FUNCTION contract_one_center_matrix(density_matrix, integral_matrix, new_to_old) RESULT(value)
2870 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: density_matrix, integral_matrix
2871 INTEGER, DIMENSION(:), INTENT(IN) :: new_to_old
2872 REAL(kind=dp) :: value
2873
2874 INTEGER :: ibas, jbas, nbas
2875
2876 nbas = SIZE(density_matrix, 1)
2877 cpassert(SIZE(density_matrix, 2) == nbas)
2878 cpassert(SIZE(new_to_old) >= nbas)
2879 cpassert(minval(new_to_old(1:nbas)) >= 1)
2880 cpassert(maxval(new_to_old(1:nbas)) <= SIZE(integral_matrix, 1))
2881 cpassert(maxval(new_to_old(1:nbas)) <= SIZE(integral_matrix, 2))
2882
2883 value = 0.0_dp
2884 DO jbas = 1, nbas
2885 DO ibas = 1, nbas
2886 value = value + density_matrix(ibas, jbas)* &
2887 integral_matrix(new_to_old(ibas), new_to_old(jbas))
2888 END DO
2889 END DO
2890
2891 END FUNCTION contract_one_center_matrix
2892
2893! **************************************************************************************************
2894!> \brief Replicate a distributed real-space field for atom-grid interpolation.
2895!> \param local_values ...
2896!> \param pw_grid ...
2897!> \param group ...
2898!> \param global_values ...
2899! **************************************************************************************************
2900 SUBROUTINE gather_native_grid_field(local_values, pw_grid, group, global_values)
2901 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN), &
2902 POINTER :: local_values
2903 TYPE(pw_grid_type), INTENT(IN), POINTER :: pw_grid
2904 TYPE(mp_para_env_type), INTENT(IN) :: group
2905 REAL(kind=dp), DIMENSION(:, :, :), INTENT(OUT), &
2906 POINTER :: global_values
2907
2908 INTEGER, DIMENSION(2, 3) :: bo
2909
2910 cpassert(ASSOCIATED(local_values))
2911 cpassert(ASSOCIATED(pw_grid))
2912 cpassert(.NOT. ASSOCIATED(global_values))
2913 bo = pw_grid%bounds_local
2914 ALLOCATE (global_values(pw_grid%bounds(1, 1):pw_grid%bounds(2, 1), &
2915 pw_grid%bounds(1, 2):pw_grid%bounds(2, 2), &
2916 pw_grid%bounds(1, 3):pw_grid%bounds(2, 3)))
2917 global_values = 0.0_dp
2918 global_values(bo(1, 1):bo(2, 1), bo(1, 2):bo(2, 2), bo(1, 3):bo(2, 3)) = &
2919 local_values
2920 CALL group%sum(global_values)
2921
2922 END SUBROUTINE gather_native_grid_field
2923
2924! **************************************************************************************************
2925!> \brief Build the tensor-product interpolation stencil for one Cartesian point.
2926!> \param stencil ...
2927!> \param pw_grid ...
2928!> \param cell ...
2929!> \param point ...
2930!> \param wrap_auxiliary_cell wrap all auxiliary-grid directions
2931!> \param indices_only skip interpolation weights when only grid indices are needed
2932! **************************************************************************************************
2933 SUBROUTINE create_native_grid_interpolation_stencil(stencil, pw_grid, cell, point, &
2934 wrap_auxiliary_cell, indices_only)
2935 TYPE(native_grid_interpolation_stencil_type), &
2936 INTENT(OUT) :: stencil
2937 TYPE(pw_grid_type), INTENT(IN), POINTER :: pw_grid
2938 TYPE(cell_type), INTENT(IN), POINTER :: cell
2939 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: point
2940 LOGICAL, INTENT(IN), OPTIONAL :: wrap_auxiliary_cell, indices_only
2941
2942 INTEGER :: idir, inode, relative_index
2943 INTEGER, DIMENSION(3) :: base
2944 LOGICAL :: build_weights, wrap_grid
2945 REAL(kind=dp), DIMENSION(3) :: fraction, relative
2946
2947 cpassert(ASSOCIATED(pw_grid))
2948 cpassert(ASSOCIATED(cell))
2949 stencil%active = .false.
2950 stencil%valid = .false.
2951 stencil%weight = 0.0_dp
2952 build_weights = .true.
2953 IF (PRESENT(indices_only)) build_weights = .NOT. indices_only
2954 wrap_grid = any(cell%perd == 0)
2955 IF (PRESENT(wrap_auxiliary_cell)) wrap_grid = wrap_auxiliary_cell
2956 relative = matmul(pw_grid%dh_inv, point)
2957 DO idir = 1, 3
2958 IF (cell%perd(idir) == 1 .OR. wrap_grid) THEN
2959 relative(idir) = modulo(relative(idir), real(pw_grid%npts(idir), kind=dp))
2960 ELSE IF (relative(idir) <= -real(native_grid_interp_offset_max, dp) .OR. &
2961 relative(idir) >= real(pw_grid%npts(idir) - &
2962 native_grid_interp_offset_min, kind=dp)) THEN
2963 RETURN
2964 END IF
2965 base(idir) = floor(relative(idir))
2966 fraction(idir) = relative(idir) - real(base(idir), kind=dp)
2967 IF (build_weights) THEN
2968 CALL native_grid_lagrange_weights(fraction(idir), stencil%weight(:, idir))
2969 END IF
2970 DO inode = 1, native_grid_interp_npts
2971 relative_index = base(idir) + native_grid_interp_offset_min + inode - 1
2972 IF (cell%perd(idir) == 1 .OR. wrap_grid) THEN
2973 relative_index = modulo(relative_index, pw_grid%npts(idir))
2974 stencil%valid(inode, idir) = .true.
2975 ELSE IF (relative_index >= 0 .AND. relative_index < pw_grid%npts(idir)) THEN
2976 stencil%valid(inode, idir) = .true.
2977 END IF
2978 stencil%relative_index(inode, idir) = relative_index
2979 END DO
2980 END DO
2981 stencil%active = .true.
2982
2983 END SUBROUTINE create_native_grid_interpolation_stencil
2984
2985! **************************************************************************************************
2986!> \brief Interpolate density, gradient, and kinetic-density fields in one stencil traversal.
2987!> \param density ...
2988!> \param grad_x ...
2989!> \param grad_y ...
2990!> \param grad_z ...
2991!> \param kin ...
2992!> \param stencil ...
2993!> \param density_value ...
2994!> \param grad_value ...
2995!> \param kin_value ...
2996! **************************************************************************************************
2997 SUBROUTINE interpolate_native_grid_fields(density, grad_x, grad_y, grad_z, kin, stencil, &
2998 density_value, grad_value, kin_value)
2999 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: density, grad_x, grad_y, grad_z, kin
3000 TYPE(native_grid_interpolation_stencil_type), &
3001 INTENT(IN) :: stencil
3002 REAL(kind=dp), INTENT(OUT) :: density_value
3003 REAL(kind=dp), DIMENSION(3), INTENT(OUT) :: grad_value
3004 REAL(kind=dp), INTENT(OUT) :: kin_value
3005
3006 INTEGER :: inode_x, inode_y, inode_z
3007 INTEGER, DIMENSION(3) :: index, lower_bound
3008 REAL(kind=dp) :: coefficient
3009
3010 cpassert(all(shape(grad_x) == shape(density)))
3011 cpassert(all(shape(grad_y) == shape(density)))
3012 cpassert(all(shape(grad_z) == shape(density)))
3013 cpassert(all(shape(kin) == shape(density)))
3014 density_value = 0.0_dp
3015 grad_value = 0.0_dp
3016 kin_value = 0.0_dp
3017 IF (.NOT. stencil%active) RETURN
3018 lower_bound = lbound(density)
3019 DO inode_z = 1, native_grid_interp_npts
3020 IF (.NOT. stencil%valid(inode_z, 3)) cycle
3021 index(3) = lower_bound(3) + stencil%relative_index(inode_z, 3)
3022 DO inode_y = 1, native_grid_interp_npts
3023 IF (.NOT. stencil%valid(inode_y, 2)) cycle
3024 index(2) = lower_bound(2) + stencil%relative_index(inode_y, 2)
3025 DO inode_x = 1, native_grid_interp_npts
3026 IF (.NOT. stencil%valid(inode_x, 1)) cycle
3027 index(1) = lower_bound(1) + stencil%relative_index(inode_x, 1)
3028 coefficient = stencil%weight(inode_x, 1)* &
3029 stencil%weight(inode_y, 2)* &
3030 stencil%weight(inode_z, 3)
3031 density_value = density_value + coefficient*density(index(1), index(2), index(3))
3032 grad_value(1) = grad_value(1) + coefficient*grad_x(index(1), index(2), index(3))
3033 grad_value(2) = grad_value(2) + coefficient*grad_y(index(1), index(2), index(3))
3034 grad_value(3) = grad_value(3) + coefficient*grad_z(index(1), index(2), index(3))
3035 kin_value = kin_value + coefficient*kin(index(1), index(2), index(3))
3036 END DO
3037 END DO
3038 END DO
3039
3040 END SUBROUTINE interpolate_native_grid_fields
3041
3042! **************************************************************************************************
3043!> \brief Group atom-grid rows by the disjoint PW tiles touched by their interpolation stencils.
3044!> \param pw_grid ...
3045!> \param cell ...
3046!> \param points atom-grid coordinates
3047!> \param wrap_auxiliary_cell wrap all auxiliary-grid directions
3048!> \param tile_count number of tiles along each PW-grid direction
3049!> \param bin_offsets CSR offsets into bin_rows
3050!> \param bin_rows atom-grid rows touching each tile
3051! **************************************************************************************************
3052 SUBROUTINE build_native_grid_adjoint_bins(pw_grid, cell, points, wrap_auxiliary_cell, &
3053 tile_count, bin_offsets, bin_rows)
3054 TYPE(pw_grid_type), INTENT(IN), POINTER :: pw_grid
3055 TYPE(cell_type), INTENT(IN), POINTER :: cell
3056 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: points
3057 LOGICAL, INTENT(IN) :: wrap_auxiliary_cell
3058 INTEGER, DIMENSION(3), INTENT(OUT) :: tile_count
3059 INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: bin_offsets, bin_rows
3060
3061 INTEGER :: candidate, ibin, idir, inode, irow, ix, &
3062 iy, iz, nbin, ntouched, position
3063 INTEGER, ALLOCATABLE, DIMENSION(:) :: bin_counts, bin_cursor, row_bin_count
3064 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: row_bins
3065 INTEGER, DIMENSION(3) :: touched_count
3066 INTEGER, DIMENSION(& native_grid_adjoint_max_tiles_per_direction, 3) :: touched_tiles
3067 TYPE(native_grid_interpolation_stencil_type) :: stencil
3068
3069 cpassert(ASSOCIATED(pw_grid))
3070 cpassert(ASSOCIATED(cell))
3071 cpassert(SIZE(points, 1) == 3)
3072 tile_count = (pw_grid%npts + native_grid_adjoint_tile_edge - 1)/ &
3073 native_grid_adjoint_tile_edge
3074 nbin = product(tile_count)
3075 ALLOCATE (row_bin_count(SIZE(points, 2)), &
3076 row_bins(native_grid_adjoint_max_bins_per_row, SIZE(points, 2)))
3077 row_bin_count = 0
3078 row_bins = 0
3079!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) &
3080!$OMP PRIVATE(candidate, idir, inode, ix, iy, iz, ntouched, stencil, &
3081!$OMP touched_count, touched_tiles) &
3082!$OMP SHARED(cell, points, pw_grid, row_bin_count, row_bins, tile_count, wrap_auxiliary_cell)
3083 DO irow = 1, SIZE(points, 2)
3084 CALL create_native_grid_interpolation_stencil( &
3085 stencil, pw_grid, cell, points(:, irow), wrap_auxiliary_cell, indices_only=.true.)
3086 IF (.NOT. stencil%active) cycle
3087 touched_count = 0
3088 touched_tiles = 0
3089 DO idir = 1, 3
3090 DO inode = 1, native_grid_interp_npts
3091 IF (.NOT. stencil%valid(inode, idir)) cycle
3092 candidate = stencil%relative_index(inode, idir)/native_grid_adjoint_tile_edge + 1
3093 IF (all(touched_tiles(1:touched_count(idir), idir) /= candidate)) THEN
3094 touched_count(idir) = touched_count(idir) + 1
3095 cpassert(touched_count(idir) <= SIZE(touched_tiles, 1))
3096 touched_tiles(touched_count(idir), idir) = candidate
3097 END IF
3098 END DO
3099 END DO
3100 ntouched = 0
3101 DO iz = 1, touched_count(3)
3102 DO iy = 1, touched_count(2)
3103 DO ix = 1, touched_count(1)
3104 ntouched = ntouched + 1
3105 cpassert(ntouched <= native_grid_adjoint_max_bins_per_row)
3106 row_bins(ntouched, irow) = 1 + touched_tiles(ix, 1) - 1 + tile_count(1)*( &
3107 touched_tiles(iy, 2) - 1 + tile_count(2)*( &
3108 touched_tiles(iz, 3) - 1))
3109 END DO
3110 END DO
3111 END DO
3112 row_bin_count(irow) = ntouched
3113 END DO
3114!$OMP END PARALLEL DO
3115
3116 ALLOCATE (bin_counts(nbin), bin_cursor(nbin), bin_offsets(nbin + 1))
3117 bin_counts = 0
3118 DO irow = 1, SIZE(points, 2)
3119 DO ibin = 1, row_bin_count(irow)
3120 bin_counts(row_bins(ibin, irow)) = bin_counts(row_bins(ibin, irow)) + 1
3121 END DO
3122 END DO
3123 bin_offsets(1) = 1
3124 DO ibin = 1, nbin
3125 bin_offsets(ibin + 1) = bin_offsets(ibin) + bin_counts(ibin)
3126 END DO
3127 ALLOCATE (bin_rows(bin_offsets(nbin + 1) - 1))
3128 bin_cursor(:) = bin_offsets(1:nbin)
3129 DO irow = 1, SIZE(points, 2)
3130 DO ibin = 1, row_bin_count(irow)
3131 candidate = row_bins(ibin, irow)
3132 position = bin_cursor(candidate)
3133 bin_rows(position) = irow
3134 bin_cursor(candidate) = position + 1
3135 END DO
3136 END DO
3137 DEALLOCATE (bin_counts, bin_cursor, row_bin_count, row_bins)
3138
3139 END SUBROUTINE build_native_grid_adjoint_bins
3140
3141! **************************************************************************************************
3142!> \brief Return the inclusive PW-grid bounds owned by one linear tile index.
3143!> \param pw_grid ...
3144!> \param tile_index linear tile index
3145!> \param tile_count number of tiles along each PW-grid direction
3146!> \param tile_lower zero-based lower grid index
3147!> \param tile_upper zero-based upper grid index
3148! **************************************************************************************************
3149 SUBROUTINE native_grid_adjoint_tile_bounds(pw_grid, tile_index, tile_count, &
3150 tile_lower, tile_upper)
3151 TYPE(pw_grid_type), INTENT(IN), POINTER :: pw_grid
3152 INTEGER, INTENT(IN) :: tile_index
3153 INTEGER, DIMENSION(3), INTENT(IN) :: tile_count
3154 INTEGER, DIMENSION(3), INTENT(OUT) :: tile_lower, tile_upper
3155
3156 INTEGER :: linear_tile
3157 INTEGER, DIMENSION(3) :: expected_tile_count, tile_coord
3158
3159 cpassert(ASSOCIATED(pw_grid))
3160 expected_tile_count = (pw_grid%npts + native_grid_adjoint_tile_edge - 1)/ &
3161 native_grid_adjoint_tile_edge
3162 cpassert(all(tile_count == expected_tile_count))
3163 cpassert(tile_index >= 1 .AND. tile_index <= product(tile_count))
3164 linear_tile = tile_index - 1
3165 tile_coord(1) = mod(linear_tile, tile_count(1))
3166 linear_tile = linear_tile/tile_count(1)
3167 tile_coord(2) = mod(linear_tile, tile_count(2))
3168 tile_coord(3) = linear_tile/tile_count(2)
3169 tile_lower = tile_coord*native_grid_adjoint_tile_edge
3170 tile_upper = min(tile_lower + native_grid_adjoint_tile_edge - 1, pw_grid%npts - 1)
3171
3172 END SUBROUTINE native_grid_adjoint_tile_bounds
3173
3174! **************************************************************************************************
3175!> \brief Apply the primitive-field interpolation transpose inside one disjoint PW tile.
3176!> \param density ...
3177!> \param grad ...
3178!> \param kin ...
3179!> \param pw_grid ...
3180!> \param stencil ...
3181!> \param density_value ...
3182!> \param grad_value ...
3183!> \param kin_value ...
3184!> \param nchannels number of spin channels accumulated by the adjoint
3185!> \param tile_lower zero-based lower grid index owned by this call
3186!> \param tile_upper zero-based upper grid index owned by this call
3187! **************************************************************************************************
3188 SUBROUTINE add_native_grid_fields_adjoint_tile( &
3189 density, grad, kin, pw_grid, stencil, density_value, grad_value, kin_value, nchannels, &
3190 tile_lower, tile_upper)
3191 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: density
3192 REAL(kind=dp), DIMENSION(:, :, :), INTENT(INOUT) :: grad
3193 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: kin
3194 TYPE(pw_grid_type), INTENT(IN), POINTER :: pw_grid
3195 TYPE(native_grid_interpolation_stencil_type), &
3196 INTENT(IN) :: stencil
3197 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: density_value
3198 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: grad_value
3199 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: kin_value
3200 INTEGER, INTENT(IN) :: nchannels
3201 INTEGER, DIMENSION(3), INTENT(IN) :: tile_lower, tile_upper
3202
3203 INTEGER :: idir, inode_x, inode_y, inode_z, ipt, &
3204 ispin
3205 INTEGER, DIMENSION(3) :: relative_index
3206 REAL(kind=dp) :: coefficient
3207
3208 cpassert(ASSOCIATED(pw_grid))
3209 cpassert(SIZE(density, 1) == product(pw_grid%npts))
3210 cpassert(SIZE(density, 2) == nchannels)
3211 cpassert(SIZE(grad, 1) == product(pw_grid%npts))
3212 cpassert(SIZE(grad, 2) == 3)
3213 cpassert(SIZE(grad, 3) == nchannels)
3214 cpassert(SIZE(kin, 1) == product(pw_grid%npts))
3215 cpassert(SIZE(kin, 2) == nchannels)
3216 cpassert(SIZE(density_value) == 2)
3217 cpassert(all(shape(grad_value) == [3, 2]))
3218 cpassert(SIZE(kin_value) == 2)
3219 cpassert(nchannels >= 1 .AND. nchannels <= 2)
3220 cpassert(all(tile_lower >= 0))
3221 cpassert(all(tile_upper >= tile_lower))
3222 cpassert(all(tile_upper < pw_grid%npts))
3223 IF (.NOT. stencil%active) RETURN
3224 DO inode_z = 1, native_grid_interp_npts
3225 IF (.NOT. stencil%valid(inode_z, 3)) cycle
3226 relative_index(3) = stencil%relative_index(inode_z, 3)
3227 IF (relative_index(3) < tile_lower(3) .OR. relative_index(3) > tile_upper(3)) cycle
3228 DO inode_y = 1, native_grid_interp_npts
3229 IF (.NOT. stencil%valid(inode_y, 2)) cycle
3230 relative_index(2) = stencil%relative_index(inode_y, 2)
3231 IF (relative_index(2) < tile_lower(2) .OR. relative_index(2) > tile_upper(2)) cycle
3232 DO inode_x = 1, native_grid_interp_npts
3233 IF (.NOT. stencil%valid(inode_x, 1)) cycle
3234 relative_index(1) = stencil%relative_index(inode_x, 1)
3235 IF (relative_index(1) < tile_lower(1) .OR. relative_index(1) > tile_upper(1)) cycle
3236 ipt = 1 + relative_index(1) + pw_grid%npts(1)*( &
3237 relative_index(2) + pw_grid%npts(2)*relative_index(3))
3238 coefficient = stencil%weight(inode_x, 1)* &
3239 stencil%weight(inode_y, 2)* &
3240 stencil%weight(inode_z, 3)
3241 DO ispin = 1, nchannels
3242 density(ipt, ispin) = density(ipt, ispin) + coefficient*density_value(ispin)
3243 DO idir = 1, 3
3244 grad(ipt, idir, ispin) = grad(ipt, idir, ispin) + &
3245 coefficient*grad_value(idir, ispin)
3246 END DO
3247 kin(ipt, ispin) = kin(ipt, ispin) + coefficient*kin_value(ispin)
3248 END DO
3249 END DO
3250 END DO
3251 END DO
3252
3253 END SUBROUTINE add_native_grid_fields_adjoint_tile
3254
3255! **************************************************************************************************
3256!> \brief Interpolate a replicated native-grid field at a Cartesian point.
3257!> \param values ...
3258!> \param pw_grid ...
3259!> \param cell ...
3260!> \param point ...
3261!> \param wrap_auxiliary_cell wrap all auxiliary-grid directions
3262!> \return ...
3263! **************************************************************************************************
3264 FUNCTION interpolate_native_grid(values, pw_grid, cell, point, wrap_auxiliary_cell) RESULT(value)
3265 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: values
3266 TYPE(pw_grid_type), INTENT(IN), POINTER :: pw_grid
3267 TYPE(cell_type), INTENT(IN), POINTER :: cell
3268 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: point
3269 LOGICAL, INTENT(IN), OPTIONAL :: wrap_auxiliary_cell
3270 REAL(kind=dp) :: value
3271
3272 INTEGER :: corner_x, corner_y, corner_z, idir
3273 INTEGER, DIMENSION(3) :: base, index, relative_index
3274 LOGICAL :: wrap_grid
3275 REAL(kind=dp) :: coefficient
3276 REAL(kind=dp), DIMENSION(3) :: fraction, relative
3277 REAL(kind=dp), &
3278 DIMENSION(native_grid_interp_npts, 3) :: weights
3279
3280 cpassert(ASSOCIATED(pw_grid))
3281 cpassert(ASSOCIATED(cell))
3282 DO idir = 1, 3
3283 cpassert(SIZE(values, idir) == pw_grid%npts(idir))
3284 END DO
3285
3286 wrap_grid = any(cell%perd == 0)
3287 IF (PRESENT(wrap_auxiliary_cell)) wrap_grid = wrap_auxiliary_cell
3288 relative = matmul(pw_grid%dh_inv, point)
3289 DO idir = 1, 3
3290 IF (cell%perd(idir) == 1 .OR. wrap_grid) THEN
3291 relative(idir) = modulo(relative(idir), real(pw_grid%npts(idir), kind=dp))
3292 ELSE IF (relative(idir) <= -real(native_grid_interp_offset_max, dp) .OR. &
3293 relative(idir) >= real(pw_grid%npts(idir) - &
3294 native_grid_interp_offset_min, kind=dp)) THEN
3295 value = 0.0_dp
3296 RETURN
3297 END IF
3298 base(idir) = floor(relative(idir))
3299 fraction(idir) = relative(idir) - real(base(idir), kind=dp)
3300 CALL native_grid_lagrange_weights(fraction(idir), weights(:, idir))
3301 END DO
3302
3303 value = 0.0_dp
3304 DO corner_z = native_grid_interp_offset_min, native_grid_interp_offset_max
3305 DO corner_y = native_grid_interp_offset_min, native_grid_interp_offset_max
3306 DO corner_x = native_grid_interp_offset_min, native_grid_interp_offset_max
3307 relative_index = base + [corner_x, corner_y, corner_z]
3308 coefficient = weights(corner_x - native_grid_interp_offset_min + 1, 1)* &
3309 weights(corner_y - native_grid_interp_offset_min + 1, 2)* &
3310 weights(corner_z - native_grid_interp_offset_min + 1, 3)
3311 DO idir = 1, 3
3312 IF (cell%perd(idir) == 1 .OR. wrap_grid) THEN
3313 relative_index(idir) = modulo(relative_index(idir), pw_grid%npts(idir))
3314 ELSE IF (relative_index(idir) < 0 .OR. &
3315 relative_index(idir) >= pw_grid%npts(idir)) THEN
3316 coefficient = 0.0_dp
3317 END IF
3318 END DO
3319 IF (coefficient == 0.0_dp) cycle
3320 index = lbound(values) + relative_index
3321 value = value + coefficient*values(index(1), index(2), index(3))
3322 END DO
3323 END DO
3324 END DO
3325
3326 END FUNCTION interpolate_native_grid
3327
3328! **************************************************************************************************
3329!> \brief Return the Cartesian gradient of native-grid interpolation at one point.
3330!> \param values ...
3331!> \param pw_grid ...
3332!> \param cell ...
3333!> \param point ...
3334!> \param wrap_auxiliary_cell wrap all auxiliary-grid directions
3335!> \return ...
3336! **************************************************************************************************
3337 FUNCTION interpolate_native_grid_gradient(values, pw_grid, cell, point, &
3338 wrap_auxiliary_cell) RESULT(gradient)
3339 REAL(kind=dp), DIMENSION(:, :, :), INTENT(IN) :: values
3340 TYPE(pw_grid_type), INTENT(IN), POINTER :: pw_grid
3341 TYPE(cell_type), INTENT(IN), POINTER :: cell
3342 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: point
3343 LOGICAL, INTENT(IN), OPTIONAL :: wrap_auxiliary_cell
3344 REAL(kind=dp), DIMENSION(3) :: gradient
3345
3346 INTEGER :: corner_x, corner_y, corner_z, idir, jdir
3347 INTEGER, DIMENSION(3) :: base, corner_index, index, relative_index
3348 LOGICAL :: wrap_grid
3349 REAL(kind=dp) :: coefficient
3350 REAL(kind=dp), DIMENSION(3) :: fraction, gradient_relative, relative
3351 REAL(kind=dp), &
3352 DIMENSION(native_grid_interp_npts, 3) :: derivative_weights, weights
3353
3354 cpassert(ASSOCIATED(pw_grid))
3355 cpassert(ASSOCIATED(cell))
3356 DO idir = 1, 3
3357 cpassert(SIZE(values, idir) == pw_grid%npts(idir))
3358 END DO
3359
3360 wrap_grid = any(cell%perd == 0)
3361 IF (PRESENT(wrap_auxiliary_cell)) wrap_grid = wrap_auxiliary_cell
3362 relative = matmul(pw_grid%dh_inv, point)
3363 DO idir = 1, 3
3364 IF (cell%perd(idir) == 1 .OR. wrap_grid) THEN
3365 relative(idir) = modulo(relative(idir), real(pw_grid%npts(idir), kind=dp))
3366 ELSE IF (relative(idir) <= -real(native_grid_interp_offset_max, dp) .OR. &
3367 relative(idir) >= real(pw_grid%npts(idir) - &
3368 native_grid_interp_offset_min, kind=dp)) THEN
3369 gradient = 0.0_dp
3370 RETURN
3371 END IF
3372 base(idir) = floor(relative(idir))
3373 fraction(idir) = relative(idir) - real(base(idir), kind=dp)
3374 CALL native_grid_lagrange_weights( &
3375 fraction(idir), weights(:, idir), derivative_weights(:, idir))
3376 END DO
3377
3378 gradient_relative = 0.0_dp
3379 DO corner_z = native_grid_interp_offset_min, native_grid_interp_offset_max
3380 DO corner_y = native_grid_interp_offset_min, native_grid_interp_offset_max
3381 DO corner_x = native_grid_interp_offset_min, native_grid_interp_offset_max
3382 relative_index = base + [corner_x, corner_y, corner_z]
3383 corner_index = [corner_x, corner_y, corner_z] - &
3384 native_grid_interp_offset_min + 1
3385 DO idir = 1, 3
3386 IF (cell%perd(idir) == 1 .OR. wrap_grid) THEN
3387 relative_index(idir) = modulo(relative_index(idir), pw_grid%npts(idir))
3388 ELSE IF (relative_index(idir) < 0 .OR. &
3389 relative_index(idir) >= pw_grid%npts(idir)) THEN
3390 EXIT
3391 END IF
3392 END DO
3393 IF (idir <= 3) cycle
3394 index = lbound(values) + relative_index
3395 DO idir = 1, 3
3396 coefficient = 1.0_dp
3397 DO jdir = 1, 3
3398 IF (jdir == idir) THEN
3399 coefficient = coefficient*derivative_weights( &
3400 corner_index(jdir), jdir)
3401 ELSE
3402 coefficient = coefficient*weights( &
3403 corner_index(jdir), jdir)
3404 END IF
3405 END DO
3406 gradient_relative(idir) = gradient_relative(idir) + &
3407 coefficient*values(index(1), index(2), index(3))
3408 END DO
3409 END DO
3410 END DO
3411 END DO
3412 gradient = matmul(transpose(pw_grid%dh_inv), gradient_relative)
3413
3414 END FUNCTION interpolate_native_grid_gradient
3415
3416! **************************************************************************************************
3417!> \brief Build the value and optional derivative weights for the native-grid interpolation.
3418!> \param fraction fractional coordinate between two grid points
3419!> \param weights interpolation weights
3420!> \param derivative_weights optional derivatives with respect to fraction
3421! **************************************************************************************************
3422 SUBROUTINE native_grid_lagrange_weights(fraction, weights, derivative_weights)
3423 REAL(kind=dp), INTENT(IN) :: fraction
3424 REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: weights
3425 REAL(kind=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: derivative_weights
3426
3427 INTEGER :: inode, jnode, knode, npoints
3428 REAL(kind=dp) :: denominator, derivative, numerator
3429 REAL(kind=dp), DIMENSION(SIZE(weights)) :: nodes
3430
3431 npoints = SIZE(weights)
3432 IF (PRESENT(derivative_weights)) THEN
3433 cpassert(SIZE(derivative_weights) == npoints)
3434 END IF
3435 DO inode = 1, npoints
3436 nodes(inode) = real(native_grid_interp_offset_min + inode - 1, dp)
3437 END DO
3438 DO inode = 1, npoints
3439 denominator = 1.0_dp
3440 numerator = 1.0_dp
3441 DO jnode = 1, npoints
3442 IF (jnode == inode) cycle
3443 denominator = denominator*(nodes(inode) - nodes(jnode))
3444 numerator = numerator*(fraction - nodes(jnode))
3445 END DO
3446 weights(inode) = numerator/denominator
3447 IF (PRESENT(derivative_weights)) THEN
3448 derivative = 0.0_dp
3449 DO knode = 1, npoints
3450 IF (knode == inode) cycle
3451 numerator = 1.0_dp
3452 DO jnode = 1, npoints
3453 IF (jnode == inode .OR. jnode == knode) cycle
3454 numerator = numerator*(fraction - nodes(jnode))
3455 END DO
3456 derivative = derivative + numerator/denominator
3457 END DO
3458 derivative_weights(inode) = derivative
3459 END IF
3460 END DO
3461
3462 END SUBROUTINE native_grid_lagrange_weights
3463
3464! **************************************************************************************************
3465!> \brief Conservatively bound periodic images that can intersect a Cartesian support sphere.
3466!> \param cell current cell and periodicity
3467!> \param fractional target-minus-source displacement in lattice coordinates
3468!> \param cutoff current, density-dependent Cartesian support radius
3469!> \param image_shell original enclosing image shell
3470!> \param lower inclusive lower image indices
3471!> \param upper inclusive upper image indices (may be smaller than lower)
3472! **************************************************************************************************
3473 SUBROUTINE atom_grid_image_bounds(cell, fractional, cutoff, image_shell, lower, upper)
3474 TYPE(cell_type), INTENT(IN), POINTER :: cell
3475 REAL(dp), INTENT(IN) :: fractional(3), cutoff
3476 INTEGER, INTENT(IN) :: image_shell(3)
3477 INTEGER, INTENT(OUT) :: lower(3), upper(3)
3478
3479 INTEGER :: center, idir
3480 REAL(dp) :: condition_bound, padding, radius
3481
3482 ! |(H^-1 d)_i| <= cutoff ||(H^-1)_i|| also holds for triclinic cells.
3483 ! Inflate the bounds for cancellation in fractional coordinates and the inverse cell.
3484 condition_bound = max(1.0_dp, sum(abs(cell%hmat))*sum(abs(cell%h_inv)))
3485 lower = 0
3486 upper = 0
3487 DO idir = 1, 3
3488 IF (cell%perd(idir) == 0) cycle
3489 center = nint(fractional(idir))
3490 radius = cutoff*sqrt(sum(cell%h_inv(idir, :)**2))
3491 padding = 128.0_dp*epsilon(1.0_dp)*condition_bound* &
3492 max(1.0_dp, maxval(abs(fractional)), radius)
3493 lower(idir) = max(center - image_shell(idir), ceiling(fractional(idir) - radius - padding))
3494 upper(idir) = min(center + image_shell(idir), floor(fractional(idir) + radius + padding))
3495 END DO
3496 END SUBROUTINE atom_grid_image_bounds
3497
3498! **************************************************************************************************
3499!> \brief Express one radial node's first and second derivatives as linear combinations of the
3500!> local node values used by the C2-continuous quintic Hermite interpolant.
3501!> \param grid_atom radial source grid
3502!> \param descending whether radial nodes are stored in descending order
3503!> \param node logical index of the node whose derivatives are required
3504!> \param logical_start first logical index represented by the coefficient arrays
3505!> \param slope coefficients of the first derivative
3506!> \param curvature coefficients of the second derivative
3507! **************************************************************************************************
3508 SUBROUTINE radial_node_derivative_coefficients( &
3509 grid_atom, descending, node, logical_start, slope, curvature)
3510 TYPE(grid_atom_type), POINTER :: grid_atom
3511 LOGICAL, INTENT(IN) :: descending
3512 INTEGER, INTENT(IN) :: node, logical_start
3513 REAL(dp), DIMENSION(4), INTENT(OUT) :: slope, curvature
3514
3515 INTEGER :: center, hi, lo, n
3516 REAL(dp) :: h_hi, h_lo, x_center, x_hi, x_lo
3517
3518 n = grid_atom%nr
3519 slope = 0.0_dp
3520 curvature = 0.0_dp
3521 IF (node == 1) THEN
3522 lo = 1
3523 hi = 2
3524 ELSE IF (node == n) THEN
3525 lo = n - 1
3526 hi = n
3527 ELSE
3528 lo = node - 1
3529 hi = node + 1
3530 END IF
3531 x_lo = grid_atom%rad(merge(n + 1 - lo, lo, descending))
3532 x_hi = grid_atom%rad(merge(n + 1 - hi, hi, descending))
3533 slope(lo - logical_start + 1) = -1.0_dp/(x_hi - x_lo)
3534 slope(hi - logical_start + 1) = 1.0_dp/(x_hi - x_lo)
3535
3536 IF (n < 3) RETURN
3537 center = min(max(node, 2), n - 1)
3538 lo = center - 1
3539 hi = center + 1
3540 x_lo = grid_atom%rad(merge(n + 1 - lo, lo, descending))
3541 x_center = grid_atom%rad(merge(n + 1 - center, center, descending))
3542 x_hi = grid_atom%rad(merge(n + 1 - hi, hi, descending))
3543 h_lo = x_center - x_lo
3544 h_hi = x_hi - x_center
3545 curvature(lo - logical_start + 1) = 2.0_dp/(h_lo*(h_lo + h_hi))
3546 curvature(center - logical_start + 1) = &
3547 -2.0_dp*(1.0_dp/h_lo + 1.0_dp/h_hi)/(h_lo + h_hi)
3548 curvature(hi - logical_start + 1) = 2.0_dp/(h_hi*(h_lo + h_hi))
3549
3550 END SUBROUTINE radial_node_derivative_coefficients
3551
3552! **************************************************************************************************
3553!> \brief Build mutually consistent value and Cartesian-derivative weights for interpolation from
3554!> a GAPW radial/Lebedev atom grid to one arbitrary local displacement.
3555!> \param grid_atom source atom grid
3556!> \param harmonics source spherical harmonics
3557!> \param displacement point minus source-image coordinate
3558!> \param cutoff compact support radius
3559!> \param radial_indices active radial nodes
3560!> \param radial_weights interpolation weights on the active radial nodes
3561!> \param radial_derivative_weights optional radial derivatives of the interpolation weights
3562!> \param nradial number of active radial nodes
3563!> \param angular_weights angular interpolation weights
3564!> \param angular_derivative_weights optional Cartesian derivatives of the angular weights
3565!> \param active whether the point lies inside the compact support
3566! **************************************************************************************************
3567 SUBROUTINE atom_grid_interpolation_weights( &
3568 grid_atom, harmonics, displacement, cutoff, radial_indices, radial_weights, &
3569 radial_derivative_weights, nradial, angular_weights, angular_derivative_weights, active)
3570 TYPE(grid_atom_type), POINTER :: grid_atom
3571 TYPE(harmonics_atom_type), POINTER :: harmonics
3572 REAL(dp), DIMENSION(3), INTENT(IN) :: displacement
3573 REAL(dp), INTENT(IN) :: cutoff
3574 INTEGER, DIMENSION(4), INTENT(OUT) :: radial_indices
3575 REAL(dp), DIMENSION(4), INTENT(OUT) :: radial_weights
3576 REAL(dp), DIMENSION(4), INTENT(OUT), OPTIONAL :: radial_derivative_weights
3577 INTEGER, INTENT(OUT) :: nradial
3578 REAL(dp), DIMENSION(:), INTENT(OUT) :: angular_weights
3579 REAL(dp), DIMENSION(:, :), INTENT(OUT), OPTIONAL :: angular_derivative_weights
3580 LOGICAL, INTENT(OUT) :: active
3581
3582 INTEGER :: ia, ic, inode, iso, l, left, left_pos, &
3583 logical_end, logical_start, lx, ly, &
3584 lz, n, right_pos, shell_index
3585 LOGICAL :: descending
3586 REAL(dp) :: dh00, dh01, dh10, dh11, dh20, dh21, h00, &
3587 h01, h10, h11, h20, h21, h_interval, &
3588 monomial, radius, solid_derivative, t, &
3589 x1, x2
3590 REAL(dp), DIMENSION(3) :: direction
3591 REAL(dp), DIMENSION(harmonics%max_s_harm) :: angular_values
3592 REAL(dp), DIMENSION(4) :: curvature_left, curvature_right, &
3593 slope_left, slope_right
3594 REAL(dp), DIMENSION(3, harmonics%max_s_harm) :: angular_value_derivatives
3595
3596 cpassert(ASSOCIATED(grid_atom))
3597 cpassert(ASSOCIATED(harmonics))
3598 cpassert(SIZE(angular_weights) == grid_atom%ng_sphere)
3599 IF (PRESENT(angular_derivative_weights)) THEN
3600 cpassert(SIZE(angular_derivative_weights, 1) == 3)
3601 cpassert(SIZE(angular_derivative_weights, 2) == grid_atom%ng_sphere)
3602 END IF
3603
3604 radial_indices = 0
3605 radial_weights = 0.0_dp
3606 IF (PRESENT(radial_derivative_weights)) radial_derivative_weights = 0.0_dp
3607 angular_weights = 0.0_dp
3608 IF (PRESENT(angular_derivative_weights)) angular_derivative_weights = 0.0_dp
3609 nradial = 0
3610 active = .false.
3611 n = grid_atom%nr
3612 descending = grid_atom%rad(1) > grid_atom%rad(n)
3613 radius = sqrt(sum(displacement**2))
3614 IF (radius > cutoff .OR. &
3615 radius > grid_atom%rad(merge(1, n, descending))) RETURN
3616 IF (radius <= 1.0e-12_dp) RETURN
3617 direction = displacement/radius
3618
3619 cpassert(n >= 2)
3620 left = n - 1
3621 IF (radius <= grid_atom%rad(merge(n, 1, descending))) THEN
3622 left = 1
3623 ELSE
3624 DO inode = 1, n - 1
3625 IF (radius <= grid_atom%rad(merge(n - inode, inode + 1, descending))) THEN
3626 left = inode
3627 EXIT
3628 END IF
3629 END DO
3630 END IF
3631
3632 logical_start = max(1, left - 1)
3633 logical_end = min(n, left + 2)
3634 nradial = logical_end - logical_start + 1
3635 DO inode = 1, nradial
3636 radial_indices(inode) = merge(n + 2 - logical_start - inode, &
3637 logical_start + inode - 1, descending)
3638 END DO
3639 left_pos = left - logical_start + 1
3640 right_pos = left_pos + 1
3641 x1 = grid_atom%rad(radial_indices(left_pos))
3642 x2 = grid_atom%rad(radial_indices(right_pos))
3643 h_interval = x2 - x1
3644 t = (radius - x1)/h_interval
3645
3646 h00 = 1.0_dp - 10.0_dp*t**3 + 15.0_dp*t**4 - 6.0_dp*t**5
3647 h10 = t - 6.0_dp*t**3 + 8.0_dp*t**4 - 3.0_dp*t**5
3648 h20 = 0.5_dp*(t**2 - 3.0_dp*t**3 + 3.0_dp*t**4 - t**5)
3649 h01 = 10.0_dp*t**3 - 15.0_dp*t**4 + 6.0_dp*t**5
3650 h11 = -4.0_dp*t**3 + 7.0_dp*t**4 - 3.0_dp*t**5
3651 h21 = 0.5_dp*(t**3 - 2.0_dp*t**4 + t**5)
3652
3653 CALL radial_node_derivative_coefficients( &
3654 grid_atom, descending, left, logical_start, slope_left, curvature_left)
3655 CALL radial_node_derivative_coefficients( &
3656 grid_atom, descending, left + 1, logical_start, slope_right, curvature_right)
3657 radial_weights(left_pos) = radial_weights(left_pos) + h00
3658 radial_weights(right_pos) = radial_weights(right_pos) + h01
3659 radial_weights(1:nradial) = radial_weights(1:nradial) + &
3660 h_interval*(h10*slope_left(1:nradial) + &
3661 h11*slope_right(1:nradial)) + &
3662 h_interval**2*(h20*curvature_left(1:nradial) + &
3663 h21*curvature_right(1:nradial))
3664 IF (PRESENT(radial_derivative_weights)) THEN
3665 dh00 = -30.0_dp*t**2 + 60.0_dp*t**3 - 30.0_dp*t**4
3666 dh10 = 1.0_dp - 18.0_dp*t**2 + 32.0_dp*t**3 - 15.0_dp*t**4
3667 dh20 = 0.5_dp*(2.0_dp*t - 9.0_dp*t**2 + 12.0_dp*t**3 - 5.0_dp*t**4)
3668 dh01 = 30.0_dp*t**2 - 60.0_dp*t**3 + 30.0_dp*t**4
3669 dh11 = -12.0_dp*t**2 + 28.0_dp*t**3 - 15.0_dp*t**4
3670 dh21 = 0.5_dp*(3.0_dp*t**2 - 8.0_dp*t**3 + 5.0_dp*t**4)
3671 radial_derivative_weights(left_pos) = &
3672 radial_derivative_weights(left_pos) + dh00/h_interval
3673 radial_derivative_weights(right_pos) = &
3674 radial_derivative_weights(right_pos) + dh01/h_interval
3675 radial_derivative_weights(1:nradial) = radial_derivative_weights(1:nradial) + &
3676 dh10*slope_left(1:nradial) + &
3677 dh11*slope_right(1:nradial) + &
3678 h_interval*(dh20*curvature_left(1:nradial) + &
3679 dh21*curvature_right(1:nradial))
3680 END IF
3681
3682 angular_values = 0.0_dp
3683 IF (PRESENT(angular_derivative_weights)) angular_value_derivatives = 0.0_dp
3684 DO iso = 1, harmonics%max_s_harm
3685 l = indso(1, iso)
3686 CALL y_lm(direction, angular_values(iso), l, indso(2, iso))
3687 IF (l == 0 .OR. .NOT. PRESENT(angular_derivative_weights)) cycle
3688 shell_index = iso - nsoset(l - 1)
3689 DO ic = 1, nco(l)
3690 lx = indco(1, ic + ncoset(l - 1))
3691 ly = indco(2, ic + ncoset(l - 1))
3692 lz = indco(3, ic + ncoset(l - 1))
3693 IF (lx > 0) THEN
3694 monomial = real(lx, dp)*direction(1)**(lx - 1)* &
3695 direction(2)**ly*direction(3)**lz
3696 angular_value_derivatives(1, iso) = angular_value_derivatives(1, iso) + &
3697 orbtramat(l)%slm(shell_index, ic)*monomial
3698 END IF
3699 IF (ly > 0) THEN
3700 monomial = direction(1)**lx*real(ly, dp)*direction(2)**(ly - 1)* &
3701 direction(3)**lz
3702 angular_value_derivatives(2, iso) = angular_value_derivatives(2, iso) + &
3703 orbtramat(l)%slm(shell_index, ic)*monomial
3704 END IF
3705 IF (lz > 0) THEN
3706 monomial = direction(1)**lx*direction(2)**ly* &
3707 REAL(lz, dp)*direction(3)**(lz - 1)
3708 angular_value_derivatives(3, iso) = angular_value_derivatives(3, iso) + &
3709 orbtramat(l)%slm(shell_index, ic)*monomial
3710 END IF
3711 END DO
3712 DO inode = 1, 3
3713 solid_derivative = angular_value_derivatives(inode, iso)
3714 angular_value_derivatives(inode, iso) = (solid_derivative - &
3715 REAL(l, dp)*angular_values(iso)* &
3716 direction(inode))/radius
3717 END DO
3718 END DO
3719
3720 DO ia = 1, grid_atom%ng_sphere
3721 angular_weights(ia) = grid_atom%wa(ia)* &
3722 dot_product(harmonics%slm(ia, 1:harmonics%max_s_harm), &
3723 angular_values)
3724 IF (PRESENT(angular_derivative_weights)) THEN
3725 DO inode = 1, 3
3726 angular_derivative_weights(inode, ia) = grid_atom%wa(ia)* &
3727 dot_product(harmonics%slm(ia, 1:harmonics%max_s_harm), &
3728 angular_value_derivatives(inode, :))
3729 END DO
3730 END IF
3731 END DO
3732 active = .true.
3733
3734 END SUBROUTINE atom_grid_interpolation_weights
3735
3736! **************************************************************************************************
3737!> \brief Interpolate hard-minus-soft primitive fields and their spatial derivatives from one
3738!> GAPW atom grid.
3739!> \param grid_atom radial and angular source grid
3740!> \param harmonics spherical-harmonic representation of the source grid
3741!> \param displacement target point relative to the source-atom image
3742!> \param cutoff compact support radius of the source fields
3743!> \param nspins number of spin channels
3744!> \param rho_h hard one-center density values
3745!> \param rho_s soft one-center density values
3746!> \param drho_h hard one-center density-gradient values
3747!> \param drho_s soft one-center density-gradient values
3748!> \param tau_h hard one-center kinetic-energy-density values
3749!> \param tau_s soft one-center kinetic-energy-density values
3750!> \param density interpolated hard-minus-soft density
3751!> \param gradient interpolated hard-minus-soft density gradient
3752!> \param kin interpolated hard-minus-soft kinetic-energy density
3753!> \param density_spatial Cartesian derivatives of density
3754!> \param gradient_spatial Cartesian derivatives of the density gradient
3755!> \param kin_spatial Cartesian derivatives of the kinetic-energy density
3756!> \param calculate_spatial evaluate spatial derivatives (default true), otherwise return zeros
3757! **************************************************************************************************
3758 SUBROUTINE interpolate_gapw_atom_grid_fields( &
3759 grid_atom, harmonics, displacement, cutoff, nspins, rho_h, rho_s, drho_h, drho_s, &
3760 tau_h, tau_s, density, gradient, kin, density_spatial, gradient_spatial, kin_spatial, calculate_spatial)
3761 TYPE(grid_atom_type), POINTER :: grid_atom
3762 TYPE(harmonics_atom_type), POINTER :: harmonics
3763 REAL(dp), DIMENSION(3), INTENT(IN) :: displacement
3764 REAL(dp), INTENT(IN) :: cutoff
3765 INTEGER, INTENT(IN) :: nspins
3766 REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: rho_h, rho_s
3767 REAL(dp), DIMENSION(:, :, :, :), INTENT(IN) :: drho_h, drho_s
3768 REAL(dp), DIMENSION(:, :, :), INTENT(IN) :: tau_h, tau_s
3769 REAL(dp), DIMENSION(2), INTENT(OUT) :: density
3770 REAL(dp), DIMENSION(3, 2), INTENT(OUT) :: gradient
3771 REAL(dp), DIMENSION(2), INTENT(OUT) :: kin
3772 REAL(dp), DIMENSION(3, 2), INTENT(OUT) :: density_spatial
3773 REAL(dp), DIMENSION(3, 3, 2), INTENT(OUT) :: gradient_spatial
3774 REAL(dp), DIMENSION(3, 2), INTENT(OUT) :: kin_spatial
3775 LOGICAL, INTENT(IN), OPTIONAL :: calculate_spatial
3776
3777 INTEGER :: ia, idir, inode, ir, ispin, jdir, nradial
3778 INTEGER, DIMENSION(4) :: radial_indices
3779 LOGICAL :: active, need_spatial
3780 REAL(dp) :: derivative_weight, weight
3781 REAL(dp), DIMENSION(grid_atom%ng_sphere) :: angular_weights
3782 REAL(dp), DIMENSION(4) :: radial_derivative_weights, radial_weights
3783 REAL(dp), DIMENSION(3, grid_atom%ng_sphere) :: angular_derivative_weights
3784
3785 density = 0.0_dp
3786 gradient = 0.0_dp
3787 kin = 0.0_dp
3788 density_spatial = 0.0_dp
3789 gradient_spatial = 0.0_dp
3790 kin_spatial = 0.0_dp
3791 need_spatial = .true.
3792 IF (PRESENT(calculate_spatial)) need_spatial = calculate_spatial
3793 IF (need_spatial) THEN
3794 CALL atom_grid_interpolation_weights( &
3795 grid_atom, harmonics, displacement, cutoff, radial_indices, radial_weights, &
3796 radial_derivative_weights, nradial, angular_weights, angular_derivative_weights, active)
3797 ELSE
3798 CALL atom_grid_interpolation_weights( &
3799 grid_atom, harmonics, displacement, cutoff, radial_indices, radial_weights, &
3800 nradial=nradial, angular_weights=angular_weights, active=active)
3801 END IF
3802 IF (active) THEN
3803 DO inode = 1, nradial
3804 ir = radial_indices(inode)
3805 DO ia = 1, grid_atom%ng_sphere
3806 weight = radial_weights(inode)*angular_weights(ia)
3807 DO ispin = 1, nspins
3808 density(ispin) = density(ispin) + &
3809 weight*(rho_h(ia, ir, ispin) - rho_s(ia, ir, ispin))
3810 kin(ispin) = kin(ispin) + &
3811 weight*(tau_h(ia, ir, ispin) - tau_s(ia, ir, ispin))
3812 DO idir = 1, 3
3813 gradient(idir, ispin) = gradient(idir, ispin) + &
3814 weight*(drho_h(idir, ia, ir, ispin) - &
3815 drho_s(idir, ia, ir, ispin))
3816 IF (.NOT. need_spatial) cycle
3817 derivative_weight = radial_derivative_weights(inode)* &
3818 displacement(idir)/sqrt(sum(displacement**2))* &
3819 angular_weights(ia) + radial_weights(inode)* &
3820 angular_derivative_weights(idir, ia)
3821 density_spatial(idir, ispin) = density_spatial(idir, ispin) + &
3822 derivative_weight*(rho_h(ia, ir, ispin) - rho_s(ia, ir, ispin))
3823 kin_spatial(idir, ispin) = kin_spatial(idir, ispin) + &
3824 derivative_weight*(tau_h(ia, ir, ispin) - tau_s(ia, ir, ispin))
3825 DO jdir = 1, 3
3826 gradient_spatial(jdir, idir, ispin) = &
3827 gradient_spatial(jdir, idir, ispin) + derivative_weight*( &
3828 drho_h(jdir, ia, ir, ispin) - drho_s(jdir, ia, ir, ispin))
3829 END DO
3830 END DO
3831 END DO
3832 END DO
3833 END DO
3834 END IF
3835 END SUBROUTINE interpolate_gapw_atom_grid_fields
3836
3837! **************************************************************************************************
3838!> \brief Apply the exact transpose of interpolate_gapw_atom_grid_fields to one-center potentials.
3839!> \param grid_atom radial and angular source grid
3840!> \param harmonics spherical-harmonic representation of the source grid
3841!> \param displacement target point relative to the source-atom image
3842!> \param cutoff compact support radius of the source fields
3843!> \param nspins number of spin channels
3844!> \param density_adjoint model derivative with respect to density
3845!> \param gradient_adjoint model derivative with respect to the density gradient
3846!> \param kin_adjoint model derivative with respect to kinetic-energy density
3847!> \param vxc_h accumulated hard one-center density potential
3848!> \param vxc_s accumulated soft one-center density potential
3849!> \param vxg_h accumulated hard one-center density-gradient potential
3850!> \param vxg_s accumulated soft one-center density-gradient potential
3851!> \param vtau_h accumulated hard one-center kinetic-energy-density potential
3852!> \param vtau_s accumulated soft one-center kinetic-energy-density potential
3853! **************************************************************************************************
3854 SUBROUTINE add_gapw_atom_grid_interpolation_adjoint( &
3855 grid_atom, harmonics, displacement, cutoff, nspins, density_adjoint, gradient_adjoint, &
3856 kin_adjoint, vxc_h, vxc_s, vxg_h, vxg_s, vtau_h, vtau_s)
3857 TYPE(grid_atom_type), POINTER :: grid_atom
3858 TYPE(harmonics_atom_type), POINTER :: harmonics
3859 REAL(dp), DIMENSION(3), INTENT(IN) :: displacement
3860 REAL(dp), INTENT(IN) :: cutoff
3861 INTEGER, INTENT(IN) :: nspins
3862 REAL(dp), DIMENSION(2), INTENT(IN) :: density_adjoint
3863 REAL(dp), DIMENSION(3, 2), INTENT(IN) :: gradient_adjoint
3864 REAL(dp), DIMENSION(2), INTENT(IN) :: kin_adjoint
3865 REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: vxc_h, vxc_s
3866 REAL(dp), DIMENSION(:, :, :, :), INTENT(INOUT) :: vxg_h, vxg_s
3867 REAL(dp), DIMENSION(:, :, :), INTENT(INOUT) :: vtau_h, vtau_s
3868
3869 INTEGER :: ia, idir, inode, ir, ispin, nradial
3870 INTEGER, DIMENSION(4) :: radial_indices
3871 LOGICAL :: active
3872 REAL(dp) :: value
3873 REAL(dp), DIMENSION(4) :: radial_weights
3874 REAL(dp), DIMENSION(grid_atom%ng_sphere) :: angular_weights
3875
3876 CALL atom_grid_interpolation_weights( &
3877 grid_atom, harmonics, displacement, cutoff, radial_indices, radial_weights, &
3878 nradial=nradial, angular_weights=angular_weights, active=active)
3879 IF (active) THEN
3880 DO inode = 1, nradial
3881 ir = radial_indices(inode)
3882 DO ia = 1, grid_atom%ng_sphere
3883 value = radial_weights(inode)*angular_weights(ia)
3884 DO ispin = 1, nspins
3885 ! CP2K applies the hard-minus-soft sign when the two one-center
3886 ! matrices are assembled, so both stored potentials carry the
3887 ! same transpose-interpolation coefficient.
3888 vxc_h(ia, ir, ispin) = vxc_h(ia, ir, ispin) + value*density_adjoint(ispin)
3889 vxc_s(ia, ir, ispin) = vxc_s(ia, ir, ispin) + value*density_adjoint(ispin)
3890 vtau_h(ia, ir, ispin) = vtau_h(ia, ir, ispin) + value*kin_adjoint(ispin)
3891 vtau_s(ia, ir, ispin) = vtau_s(ia, ir, ispin) + value*kin_adjoint(ispin)
3892 DO idir = 1, 3
3893 vxg_h(idir, ia, ir, ispin) = vxg_h(idir, ia, ir, ispin) + &
3894 value*gradient_adjoint(idir, ispin)
3895 vxg_s(idir, ia, ir, ispin) = vxg_s(idir, ia, ir, ispin) + &
3896 value*gradient_adjoint(idir, ispin)
3897 END DO
3898 END DO
3899 END DO
3900 END DO
3901 END IF
3902 END SUBROUTINE add_gapw_atom_grid_interpolation_adjoint
3903
3904! **************************************************************************************************
3905!> \brief ...
3906!> \param qs_env ...
3907!> \param exc1 the on-body ex energy contribution
3908!> \param gradient_atom_set ...
3909! **************************************************************************************************
3910 SUBROUTINE calculate_vxc_atom_epr(qs_env, exc1, gradient_atom_set)
3911
3912 TYPE(qs_environment_type), POINTER :: qs_env
3913 REAL(dp), INTENT(INOUT) :: exc1
3914 TYPE(nablavks_atom_type), DIMENSION(:), POINTER :: gradient_atom_set
3915
3916 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_vxc_atom_epr'
3917
3918 INTEGER :: bo(2), handle, ia, iat, iatom, idir, &
3919 ikind, ir, ispin, myfun, na, natom, &
3920 nr, nspins, num_pe, on
3921 INTEGER, DIMENSION(2, 3) :: bounds
3922 INTEGER, DIMENSION(:), POINTER :: atom_list
3923 LOGICAL :: accint, donlcc, gradient_f, lsd, nlcc, &
3924 paw_atom, tau_f
3925 REAL(dp) :: agr, alpha, density_cut, exc_h, exc_s, &
3926 gradient_cut, tau_cut
3927 REAL(dp), ALLOCATABLE, DIMENSION(:) :: fw
3928 REAL(dp), DIMENSION(1, 1, 1) :: tau_d
3929 REAL(dp), DIMENSION(1, 1, 1, 1) :: rho_d
3930 REAL(dp), DIMENSION(:, :), POINTER :: rho_nlcc, weight_h, weight_s
3931 REAL(dp), DIMENSION(:, :, :), POINTER :: rho_h, rho_s, tau_h, tau_s, vtau_h, &
3932 vtau_s, vxc_h, vxc_s
3933 REAL(dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg_h, vxg_s
3934 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
3935 TYPE(dft_control_type), POINTER :: dft_control
3936 TYPE(grid_atom_type), POINTER :: grid_atom
3937 TYPE(gto_basis_set_type), POINTER :: basis_1c
3938 TYPE(harmonics_atom_type), POINTER :: harmonics
3939 TYPE(mp_para_env_type), POINTER :: para_env
3940 TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
3941 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
3942 TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
3943 TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom_set
3944 TYPE(rho_atom_type), POINTER :: rho_atom
3945 TYPE(section_vals_type), POINTER :: input, my_xc_section, xc_fun_section
3946 TYPE(tau_basis_cache_type) :: tau_basis_cache
3947 TYPE(xc_derivative_set_type) :: deriv_set
3948 TYPE(xc_rho_cflags_type) :: needs
3949 TYPE(xc_rho_set_type) :: rho_set_h, rho_set_s
3950
3951! -------------------------------------------------------------------------
3952
3953 CALL timeset(routinen, handle)
3954
3955 NULLIFY (atom_list)
3956 NULLIFY (my_kind_set)
3957 NULLIFY (atomic_kind_set)
3958 NULLIFY (grid_atom)
3959 NULLIFY (harmonics)
3960 NULLIFY (input)
3961 NULLIFY (para_env)
3962 NULLIFY (rho_atom)
3963 NULLIFY (my_rho_atom_set)
3964 NULLIFY (rho_nlcc)
3965
3966 CALL get_qs_env(qs_env=qs_env, &
3967 dft_control=dft_control, &
3968 para_env=para_env, &
3969 atomic_kind_set=atomic_kind_set, &
3970 qs_kind_set=my_kind_set, &
3971 input=input, &
3972 rho_atom_set=my_rho_atom_set)
3973
3974 nlcc = has_nlcc(my_kind_set)
3975 accint = dft_control%qs_control%gapw_control%accurate_xcint
3976
3977 my_xc_section => section_vals_get_subs_vals(input, &
3978 "PROPERTIES%LINRES%EPR%PRINT%G_TENSOR%XC")
3979 xc_fun_section => section_vals_get_subs_vals(my_xc_section, "XC_FUNCTIONAL")
3980 CALL section_vals_val_get(xc_fun_section, "_SECTION_PARAMETERS_", &
3981 i_val=myfun)
3982
3983 IF (myfun == xc_none) THEN
3984 exc1 = 0.0_dp
3985 my_rho_atom_set(:)%exc_h = 0.0_dp
3986 my_rho_atom_set(:)%exc_s = 0.0_dp
3987 ELSE
3988 CALL section_vals_val_get(my_xc_section, "DENSITY_CUTOFF", &
3989 r_val=density_cut)
3990 CALL section_vals_val_get(my_xc_section, "GRADIENT_CUTOFF", &
3991 r_val=gradient_cut)
3992 CALL section_vals_val_get(my_xc_section, "TAU_CUTOFF", &
3993 r_val=tau_cut)
3994
3995 lsd = dft_control%lsd
3996 nspins = dft_control%nspins
3997 needs = xc_functionals_get_needs(xc_fun_section, &
3998 lsd=lsd, &
3999 calc_potential=.true.)
4000
4001 ! whatever the xc, if epr_xc, drho_spin is needed
4002 needs%drho_spin = .true.
4003
4004 gradient_f = (needs%drho .OR. needs%drho_spin)
4005 tau_f = (needs%tau .OR. needs%tau_spin)
4006
4007 ! Initialize energy contribution from the one center XC terms to zero
4008 exc1 = 0.0_dp
4009
4010 ! Nullify some pointers for work-arrays
4011 NULLIFY (rho_h, drho_h, rho_s, drho_s, weight_h, weight_s)
4012 NULLIFY (vxc_h, vxc_s, vxg_h, vxg_s)
4013 NULLIFY (tau_h, tau_s)
4014 NULLIFY (vtau_h, vtau_s)
4015
4016 ! Here starts the loop over all the atoms
4017
4018 DO ikind = 1, SIZE(atomic_kind_set)
4019 CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
4020 CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom, &
4021 harmonics=harmonics, grid_atom=grid_atom)
4022 CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_1c, basis_type="GAPW_1C")
4023
4024 IF (.NOT. paw_atom) cycle
4025
4026 nr = grid_atom%nr
4027 na = grid_atom%ng_sphere
4028
4029 ! Prepare the structures needed to calculate and store the xc derivatives
4030
4031 ! Array dimension: here anly one dimensional arrays are used,
4032 ! i.e. only the first column of deriv_data is read.
4033 ! The other to dimensions are set to size equal 1
4034 bounds(1:2, 1:3) = 1
4035 bounds(2, 1) = na
4036 bounds(2, 2) = nr
4037
4038 ! set integration weights
4039 weight_h => grid_atom%weight
4040 weight_s => grid_atom%weight
4041 IF (accint) THEN
4042 on = dft_control%qs_control%gapw_control%oweights
4043 alpha = dft_control%qs_control%gapw_control%aw(ikind)
4044 IF (ASSOCIATED(grid_atom%gapw_weight_s)) THEN
4045 IF (grid_atom%gapw_weight_alpha /= alpha) DEALLOCATE (grid_atom%gapw_weight_s)
4046 END IF
4047 IF (.NOT. ASSOCIATED(grid_atom%gapw_weight_s)) THEN
4048 ALLOCATE (grid_atom%gapw_weight_s(na, nr))
4049 ALLOCATE (fw(nr))
4050 CALL calc_weight_function(fw, grid_atom%rad2, on, alpha)
4051 DO ir = 1, nr
4052 agr = 1.0_dp - fw(ir)
4053 grid_atom%gapw_weight_s(:, ir) = agr*grid_atom%weight(:, ir)
4054 END DO
4055 DEALLOCATE (fw)
4056 grid_atom%gapw_weight_alpha = alpha
4057 END IF
4058 weight_s => grid_atom%gapw_weight_s
4059 END IF
4060
4061 ! create a place where to put the derivatives
4062 CALL xc_dset_create(deriv_set, local_bounds=bounds)
4063 ! create the place where to store the argument for the functionals
4064 CALL xc_rho_set_create(rho_set_h, bounds, rho_cutoff=density_cut, &
4065 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
4066 CALL xc_rho_set_create(rho_set_s, bounds, rho_cutoff=density_cut, &
4067 drho_cutoff=gradient_cut, tau_cutoff=tau_cut)
4068
4069 ! allocate the required 3d arrays where to store rho and drho
4070 CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
4071 CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
4072
4073 CALL reallocate(rho_h, 1, na, 1, nr, 1, nspins)
4074 CALL reallocate(rho_s, 1, na, 1, nr, 1, nspins)
4075 CALL reallocate(vxc_h, 1, na, 1, nr, 1, nspins)
4076 CALL reallocate(vxc_s, 1, na, 1, nr, 1, nspins)
4077 !
4078 IF (gradient_f) THEN
4079 CALL reallocate(drho_h, 1, 4, 1, na, 1, nr, 1, nspins)
4080 CALL reallocate(drho_s, 1, 4, 1, na, 1, nr, 1, nspins)
4081 CALL reallocate(vxg_h, 1, 3, 1, na, 1, nr, 1, nspins)
4082 CALL reallocate(vxg_s, 1, 3, 1, na, 1, nr, 1, nspins)
4083 END IF
4084
4085 IF (tau_f) THEN
4086 CALL create_tau_basis_cache(tau_basis_cache, grid_atom, basis_1c, harmonics)
4087 CALL reallocate(tau_h, 1, na, 1, nr, 1, nspins)
4088 CALL reallocate(tau_s, 1, na, 1, nr, 1, nspins)
4089 CALL reallocate(vtau_h, 1, na, 1, nr, 1, nspins)
4090 CALL reallocate(vtau_s, 1, na, 1, nr, 1, nspins)
4091 END IF
4092
4093 ! NLCC: prepare rho and drho of the core charge for this KIND
4094 donlcc = .false.
4095 IF (nlcc) THEN
4096 NULLIFY (rho_nlcc)
4097 rho_nlcc => my_kind_set(ikind)%nlcc_pot
4098 IF (ASSOCIATED(rho_nlcc)) donlcc = .true.
4099 END IF
4100
4101 ! Distribute the atoms of this kind
4102
4103 num_pe = para_env%num_pe
4104 bo = get_limit(natom, para_env%num_pe, para_env%mepos)
4105
4106 DO iat = bo(1), bo(2)
4107 iatom = atom_list(iat)
4108
4109 my_rho_atom_set(iatom)%exc_h = 0.0_dp
4110 my_rho_atom_set(iatom)%exc_s = 0.0_dp
4111
4112 rho_atom => my_rho_atom_set(iatom)
4113 rho_h = 0.0_dp
4114 rho_s = 0.0_dp
4115 IF (gradient_f) THEN
4116 NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
4117 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
4118 rho_rad_s=r_s, drho_rad_h=dr_h, &
4119 drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
4120 rho_rad_s_d=r_s_d)
4121 drho_h = 0.0_dp
4122 drho_s = 0.0_dp
4123 ELSE
4124 NULLIFY (r_h, r_s)
4125 CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
4126 rho_d = 0.0_dp
4127 END IF
4128 IF (tau_f) THEN
4129 !compute tau on the grid all at once
4130 CALL calc_tau_atom(tau_h, tau_s, rho_atom, tau_basis_cache, nspins)
4131 ELSE
4132 tau_d = 0.0_dp
4133 END IF
4134
4135 DO ir = 1, nr
4136 CALL calc_rho_angular(grid_atom, harmonics, nspins, gradient_f, &
4137 ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, &
4138 r_h_d, r_s_d, drho_h, drho_s)
4139 IF (donlcc) THEN
4140 CALL calc_rho_nlcc(grid_atom, nspins, gradient_f, &
4141 ir, rho_nlcc(:, 1), rho_h, rho_s, rho_nlcc(:, 2), drho_h, drho_s)
4142 END IF
4143 END DO
4144 DO ir = 1, nr
4145 IF (tau_f) THEN
4146 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_h, na, ir)
4147 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_s, na, ir)
4148 ELSE IF (gradient_f) THEN
4149 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau_d, na, ir)
4150 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau_d, na, ir)
4151 ELSE
4152 CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, rho_d, tau_d, na, ir)
4153 CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, rho_d, tau_d, na, ir)
4154 END IF
4155 END DO
4156
4157 !-------------------!
4158 ! hard atom density !
4159 !-------------------!
4160 CALL xc_dset_zero_all(deriv_set)
4161 CALL vxc_of_r_epr(xc_fun_section, rho_set_h, deriv_set, needs, weight_h, &
4162 lsd, na, nr, exc_h, vxc_h, vxg_h, vtau_h)
4163 rho_atom%exc_h = rho_atom%exc_h + exc_h
4164
4165 !-------------------!
4166 ! soft atom density !
4167 !-------------------!
4168 CALL xc_dset_zero_all(deriv_set)
4169 CALL vxc_of_r_epr(xc_fun_section, rho_set_s, deriv_set, needs, weight_s, &
4170 lsd, na, nr, exc_s, vxc_s, vxg_s, vtau_s)
4171 rho_atom%exc_s = rho_atom%exc_s + exc_s
4172
4173 DO ispin = 1, nspins
4174 DO idir = 1, 3
4175 DO ir = 1, nr
4176 DO ia = 1, na
4177 gradient_atom_set(iatom)%nablavks_vec_rad_h(idir, ispin)%r_coef(ir, ia) = &
4178 gradient_atom_set(iatom)%nablavks_vec_rad_h(idir, ispin)%r_coef(ir, ia) &
4179 + vxg_h(idir, ia, ir, ispin)
4180 gradient_atom_set(iatom)%nablavks_vec_rad_s(idir, ispin)%r_coef(ir, ia) = &
4181 gradient_atom_set(iatom)%nablavks_vec_rad_s(idir, ispin)%r_coef(ir, ia) &
4182 + vxg_s(idir, ia, ir, ispin)
4183 END DO ! ia
4184 END DO ! ir
4185 END DO ! idir
4186 END DO ! ispin
4187
4188 ! Add contributions to the exc energy
4189
4190 exc1 = exc1 + rho_atom%exc_h - rho_atom%exc_s
4191
4192 ! Integration to get the matrix elements relative to the vxc_atom
4193 ! here the products with the primitives is done: gaVxcgb
4194 ! internal transformation to get the integral in cartesian Gaussians
4195
4196 NULLIFY (int_hh, int_ss)
4197 CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
4198 IF (gradient_f) THEN
4199 CALL gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, &
4200 grid_atom, basis_1c, harmonics, nspins)
4201 ELSE
4202 CALL gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, &
4203 grid_atom, basis_1c, harmonics, nspins)
4204 END IF
4205 IF (tau_f) THEN
4206 CALL dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, &
4207 tau_basis_cache, nspins)
4208 END IF
4209 NULLIFY (r_h, r_s, dr_h, dr_s)
4210 END DO ! iat
4211
4212 IF (tau_f) CALL release_tau_basis_cache(tau_basis_cache)
4213
4214 ! Release the xc structure used to store the xc derivatives
4215 CALL xc_dset_release(deriv_set)
4216 CALL xc_rho_set_release(rho_set_h)
4217 CALL xc_rho_set_release(rho_set_s)
4218 END DO ! ikind
4219
4220 CALL para_env%sum(exc1)
4221
4222 IF (ASSOCIATED(rho_h)) DEALLOCATE (rho_h)
4223 IF (ASSOCIATED(rho_s)) DEALLOCATE (rho_s)
4224 IF (ASSOCIATED(vxc_h)) DEALLOCATE (vxc_h)
4225 IF (ASSOCIATED(vxc_s)) DEALLOCATE (vxc_s)
4226
4227 IF (gradient_f) THEN
4228 IF (ASSOCIATED(drho_h)) DEALLOCATE (drho_h)
4229 IF (ASSOCIATED(drho_s)) DEALLOCATE (drho_s)
4230 IF (ASSOCIATED(vxg_h)) DEALLOCATE (vxg_h)
4231 IF (ASSOCIATED(vxg_s)) DEALLOCATE (vxg_s)
4232 END IF
4233
4234 IF (tau_f) THEN
4235 IF (ASSOCIATED(tau_h)) DEALLOCATE (tau_h)
4236 IF (ASSOCIATED(tau_s)) DEALLOCATE (tau_s)
4237 IF (ASSOCIATED(vtau_h)) DEALLOCATE (vtau_h)
4238 IF (ASSOCIATED(vtau_s)) DEALLOCATE (vtau_s)
4239 END IF
4240
4241 END IF !xc_none
4242
4243 CALL timestop(handle)
4244
4245 END SUBROUTINE calculate_vxc_atom_epr
4246
4247END MODULE qs_vxc_atom
4248
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
memory_utilities::reallocate
Definition memory_utilities.F:27
pw_methods::pw_axpy
Definition pw_methods.F:165
spherical_harmonics::y_lm
Definition spherical_harmonics.F:45
atom
Definition atom.F:9
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_array_utils
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Definition cp_array_utils.F:18
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::cdft_beta_constraint
integer, parameter, public cdft_beta_constraint
Definition input_constants.F:820
input_constants::cdft_magnetization_constraint
integer, parameter, public cdft_magnetization_constraint
Definition input_constants.F:820
input_constants::cdft_charge_constraint
integer, parameter, public cdft_charge_constraint
Definition input_constants.F:820
input_constants::cdft_alpha_constraint
integer, parameter, public cdft_alpha_constraint
Definition input_constants.F:820
input_constants::xc_none
integer, parameter, public xc_none
Definition input_constants.F:593
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_get_rval
real(kind=dp) function, public section_get_rval(section_vals, keyword_name)
...
Definition input_section_types.F:953
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::int_8
integer, parameter, public int_8
Definition kinds.F:54
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
lebedev
Generation of the spherical Lebedev grids. All Lebedev grids were generated with a precision of at le...
Definition lebedev.F:57
lebedev::deallocate_lebedev_grids
subroutine, public deallocate_lebedev_grids()
...
Definition lebedev.F:324
lebedev::lebedev_grid
type(oh_grid), dimension(nlg), target, public lebedev_grid
Definition lebedev.F:85
lebedev::get_number_of_lebedev_grid
integer function, public get_number_of_lebedev_grid(l, n)
Get the number of the Lebedev grid, which has the requested angular momentum quantnum number l or siz...
Definition lebedev.F:114
lebedev::init_lebedev_grids
subroutine, public init_lebedev_grids()
Load the coordinates and weights of the nonredundant Lebedev grid points.
Definition lebedev.F:344
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition orbital_pointers.F:43
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition orbital_transformation_matrices.F:47
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_env_types::pw_env_get
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition pw_env_types.F:113
pw_grid_types
Definition pw_grid_types.F:13
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_cdft_grid
Pointwise CDFT partition functions for nonuniform integration grids.
Definition qs_cdft_grid.F:11
qs_cdft_grid::cdft_point_context_create
subroutine, public cdft_point_context_create(qs_env, context, calculate_derivatives)
Initialize reusable data for pointwise CDFT partition evaluation.
Definition qs_cdft_grid.F:72
qs_cdft_grid::cdft_point_context_release
subroutine, public cdft_point_context_release(context)
Release a pointwise CDFT partition context.
Definition qs_cdft_grid.F:290
qs_cdft_grid::cdft_point_weights
subroutine, public cdft_point_weights(context, point, weights, point_derivative, atom_derivative, atomic_weights)
Evaluate CDFT weights and coordinate derivatives at one point.
Definition qs_cdft_grid.F:330
qs_cdft_types
Defines CDFT control structures.
Definition qs_cdft_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:545
qs_force_types
Definition qs_force_types.F:13
qs_grid_atom
Definition qs_grid_atom.F:8
qs_grid_atom::allocate_grid_atom
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
Definition qs_grid_atom.F:74
qs_grid_atom::create_grid_atom
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Definition qs_grid_atom.F:189
qs_harmonics_atom
Definition qs_harmonics_atom.F:8
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::has_nlcc
logical function, public has_nlcc(qs_kind_set)
finds if a given qs run needs to use nlcc
Definition qs_kind_types.F:4291
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_linres_types
Type definitiona for linear response calculations.
Definition qs_linres_types.F:12
qs_native_grid_cache
Bounded geometry-only cache for native atom-grid interpolation.
Definition qs_native_grid_cache.F:8
qs_native_grid_cache::fetch_native_grid_stencil
logical function, public fetch_native_grid_stencil(cache, row, point, stencil)
Read a stencil only on an exact coordinate match. Safe for concurrent readers.
Definition qs_native_grid_cache.F:115
qs_native_grid_cache::store_native_grid_stencil
subroutine, public store_native_grid_stencil(cache, row, point, stencil)
Store from the unique forward owner of a row, never from concurrent adjoint readers.
Definition qs_native_grid_cache.F:139
qs_native_grid_cache::native_grid_interp_npts
integer, parameter, public native_grid_interp_npts
Definition qs_native_grid_cache.F:18
qs_native_grid_cache::native_grid_interp_offset_max
integer, parameter, public native_grid_interp_offset_max
Definition qs_native_grid_cache.F:17
qs_native_grid_cache::prepare_native_grid_cache
subroutine, public prepare_native_grid_cache(cache, grid, cell, nrows, wrap, max_bytes)
Prepare outside parallel regions. Cache a bounded prefix and compute the rest normally.
Definition qs_native_grid_cache.F:56
qs_native_grid_cache::native_grid_interp_offset_min
integer, parameter, public native_grid_interp_offset_min
Definition qs_native_grid_cache.F:17
qs_rho_atom_methods
Definition qs_rho_atom_methods.F:7
qs_rho_atom_methods::replicate_rho_atom_radial
subroutine, public replicate_rho_atom_radial(para_env, rho_atom_set, qs_kind, atom_list, natom, nspins)
Replicate the radial hard/soft density data needed to evaluate one-center tails on rank-local target ...
Definition qs_rho_atom_methods.F:449
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_atom_types::get_rho_atom
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_vlocal_gb_h, ga_vlocal_gb_s, int_scr_h, int_scr_s)
...
Definition qs_rho_atom_types.F:195
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_vxc_atom_utils
Support routines for integrals of the Vxc/Fxc/Gxc potentials calculated for the atomic density in the...
Definition qs_vxc_atom_utils.F:12
qs_vxc_atom_utils::calc_rho_angular
subroutine, public calc_rho_angular(grid_atom, harmonics, nspins, grad_func, ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
...
Definition qs_vxc_atom_utils.F:219
qs_vxc_atom_utils::release_tau_basis_cache
subroutine, public release_tau_basis_cache(tau_cache)
Release precomputed GAPW meta-GGA tau factors.
Definition qs_vxc_atom_utils.F:402
qs_vxc_atom_utils::create_tau_basis_cache
subroutine, public create_tau_basis_cache(tau_cache, grid_atom, basis_1c, harmonics)
Precompute radial and angular factors for GAPW meta-GGA tau contractions.
Definition qs_vxc_atom_utils.F:330
qs_vxc_atom_utils::gavxcgb_nogc
subroutine, public gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition qs_vxc_atom_utils.F:595
qs_vxc_atom_utils::gavxcgb_gc
subroutine, public gavxcgb_gc(vxc_h, vxc_s, vxg_h, vxg_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition qs_vxc_atom_utils.F:759
qs_vxc_atom_utils::calc_weight_function
subroutine, public calc_weight_function(fun, r2, order, alpha)
Calculates the radial weight function for accurate XC integration.
Definition qs_vxc_atom_utils.F:1108
qs_vxc_atom_utils::calc_rho_nlcc
subroutine, public calc_rho_nlcc(grid_atom, nspins, grad_func, ir, rho_nlcc, rho_h, rho_s, drho_nlcc, drho_h, drho_s)
...
Definition qs_vxc_atom_utils.F:516
qs_vxc_atom_utils::dgavtaudgb
subroutine, public dgavtaudgb(vtau_h, vtau_s, int_hh, int_ss, tau_cache, nspins)
Integrates 0.5 * grad_ga .dot. (V_tau * grad_gb) on the atomic grid for meta-GGA.
Definition qs_vxc_atom_utils.F:1012
qs_vxc_atom_utils::calc_tau_atom
subroutine, public calc_tau_atom(tau_h, tau_s, rho_atom, tau_cache, nspins)
Computes tau hard and soft on the atomic grids for meta-GGA calculations.
Definition qs_vxc_atom_utils.F:430
qs_vxc_atom_utils::gapw_atom_grid_support_radius
real(dp) function, public gapw_atom_grid_support_radius(grid_atom, rho_h, rho_s, drho_h, drho_s, tau_h, tau_s)
Compact support radius of one hard-minus-soft atom-grid field.
Definition qs_vxc_atom_utils.F:158
qs_vxc_atom_utils::evaluate_nlcc_primitive_fields
subroutine, public evaluate_nlcc_primitive_fields(point, center, gth_potential, sgp_potential, rho, gradient, hessian)
Evaluate an NLCC density and its first two Cartesian derivatives.
Definition qs_vxc_atom_utils.F:70
qs_vxc_atom
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition qs_vxc_atom.F:12
qs_vxc_atom::gapw_cdft_one_center
subroutine, public gapw_cdft_one_center(qs_env, energy_only, calculate_forces, values, electronic_charge, operator_group, rho_atom_operator_set)
Add the GAPW one-center correction to CDFT values and operators.
Definition qs_vxc_atom.F:2638
qs_vxc_atom::calculate_vxc_atom
subroutine, public calculate_vxc_atom(qs_env, energy_only, exc1, adiabatic_rescale_factor, kind_set_external, rho_atom_set_external, xc_section_external, calculate_forces, composite_vxc_rho, composite_vxc_tau, composite_reference_active, direct_valence_atom_grid, atom_composite_grid)
...
Definition qs_vxc_atom.F:231
qs_vxc_atom::calculate_vxc_atom_epr
subroutine, public calculate_vxc_atom_epr(qs_env, exc1, gradient_atom_set)
...
Definition qs_vxc_atom.F:3912
skala_gpw_features
Build SKALA TorchScript feature dictionaries from CP2K GPW real-space grids.
Definition skala_gpw_features.F:11
skala_gpw_features::periodic_atom_image_partition_from_layout
subroutine, public periodic_atom_image_partition_from_layout(grid_point, image_coords, target_image, weight)
Return an image-complete periodic atom weight for a prebuilt image layout.
Definition skala_gpw_features.F:2769
skala_gpw_features::build_periodic_atom_image_layout
subroutine, public build_periodic_atom_image_layout(atom_coords, cell, target_atom, image_periodicity, image_coords, target_image)
Build the image coordinates shared by all points of one target-atom block.
Definition skala_gpw_features.F:2697
skala_gpw_features::smooth_partition_atomic_weight_scale_derivative
pure real(kind=dp) function, public smooth_partition_atomic_weight_scale_derivative(weight)
Derivative of the sparse atom-row quadrature taper with respect to partition weight.
Definition skala_gpw_features.F:3240
skala_gpw_features::smooth_atom_partition
subroutine, public smooth_atom_partition(grid_point, atom_coords, cell, weights, partition_atom_coords, distances, pair_distances)
Build Becke-like smooth atom weights for one native-grid point.
Definition skala_gpw_features.F:2432
skala_gpw_features::periodic_atom_image_partition
subroutine, public periodic_atom_image_partition(grid_point, atom_coords, cell, target_atom, weight, dweight_datom, dweight_dstrain, image_periodicity)
Return the smooth weight of one reference-cell atom in an image-complete periodic atom partition....
Definition skala_gpw_features.F:2838
skala_gpw_features::smooth_partition_atomic_weight_scale
pure real(kind=dp) function, public smooth_partition_atomic_weight_scale(weight)
Smoothly suppress a sparse atom row's internal quadrature weight at the layout cutoff.
Definition skala_gpw_features.F:3217
skala_gpw_features::skala_gpw_smooth_partition_derivatives
subroutine, public skala_gpw_smooth_partition_derivatives(grid_point, atom_coords, cell, weights, included, dweights_datom, dweights_dstrain)
Build smooth atom weights and their atom/cell deformation derivatives.
Definition skala_gpw_features.F:3044
skala_gpw_functional
Experimental CP2K-native GPW real-space-grid path for SKALA TorchScript models.
Definition skala_gpw_functional.F:11
skala_gpw_functional::skala_gapw_density_partition_soft_only
integer, parameter, public skala_gapw_density_partition_soft_only
Definition skala_gpw_functional.F:78
skala_gpw_functional::skala_gapw_atom_composite_energy
subroutine, public skala_gapw_atom_composite_energy(xc_section, group, density, grad, kin, grid_coords, grid_weights, atomic_grid_weights, atomic_grid_sizes, atomic_coords, exc, density_grad_out, grad_grad_out, kin_grad_out, grid_coord_grad_out, grid_weight_grad_out, atomic_grid_weight_grad_out, atom_coord_grad_out)
Evaluate a rank-local set of complete atom blocks and sum their SKALA energies.
Definition skala_gpw_functional.F:1157
skala_gpw_functional::build_vxc_from_feature_grads
subroutine, public build_vxc_from_feature_grads(vxc_rho, vxc_tau, rho_r, pw_pool, density_grad, grad_grad, kin_grad, xc_deriv_method_id, global_grid_layout)
Fill CP2K VXC real-space arrays from Torch feature gradients.
Definition skala_gpw_functional.F:2717
skala_gpw_functional::skala_gapw_density_partition_none
integer, parameter, public skala_gapw_density_partition_none
Definition skala_gpw_functional.F:78
skala_gpw_functional::xc_section_uses_gauxc_model
logical function, public xc_section_uses_gauxc_model(xc_section)
Return true if the GAUXC subsection requests a model evaluation.
Definition skala_gpw_functional.F:280
skala_gpw_functional::xc_section_uses_native_skala_evaluator
logical function, public xc_section_uses_native_skala_evaluator(xc_section)
Return true when SKALA must be evaluated by the CP2K-native grid machinery.
Definition skala_gpw_functional.F:165
skala_gpw_functional::skala_gapw_density_partition_hard_minus_soft
integer, parameter, public skala_gapw_density_partition_hard_minus_soft
Definition skala_gpw_functional.F:78
skala_gpw_functional::native_skala_gapw_density_partition
integer function, public native_skala_gapw_density_partition(xc_section)
Return the hard/soft GAPW one-center density partition for native SKALA.
Definition skala_gpw_functional.F:307
skala_gpw_functional::get_gauxc_section
type(section_vals_type) function, pointer, public get_gauxc_section(xc_section)
Return the first GAUXC functional subsection, if present.
Definition skala_gpw_functional.F:3037
skala_gpw_functional::skala_gapw_atom_vxc_of_r
subroutine, public skala_gapw_atom_vxc_of_r(xc_section, grid_atom, group, atom_coord, rho, drho, tau, weights, lsd, nspins, na, nr, exc, vxc, vxg, vtau, energy_only, atom_force, atom_virial)
Evaluate SKALA on a GAPW one-center atomic grid.
Definition skala_gpw_functional.F:895
skala_gpw_functional::skala_gapw_density_partition_hard_only
integer, parameter, public skala_gapw_density_partition_hard_only
Definition skala_gpw_functional.F:78
spherical_harmonics
Calculate spherical harmonics.
Definition spherical_harmonics.F:27
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
virial_types
Definition virial_types.F:13
xc_atom
Definition xc_atom.F:9
xc_atom::vxc_of_r_epr
subroutine, public vxc_of_r_epr(xc_fun_section, rho_set, deriv_set, needs, w, lsd, na, nr, exc, vxc, vxg, vtau)
Specific EPR version of vxc_of_r_new.
Definition xc_atom.F:313
xc_atom::vxc_of_r_new
subroutine, public vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, lsd, na, nr, exc, vxc, vxg, vtau, energy_only, adiabatic_rescale_factor)
...
Definition xc_atom.F:64
xc_atom::xc_rho_set_atom_update
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition xc_atom.F:523
xc_atom::fill_rho_set
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition xc_atom.F:683
xc_derivative_set_types
represent a group ofunctional derivatives
Definition xc_derivative_set_types.F:14
xc_derivative_set_types::xc_dset_zero_all
subroutine, public xc_dset_zero_all(deriv_set)
...
Definition xc_derivative_set_types.F:166
xc_derivative_set_types::xc_dset_release
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
Definition xc_derivative_set_types.F:143
xc_derivative_set_types::xc_dset_create
subroutine, public xc_dset_create(derivative_set, pw_pool, local_bounds)
creates a derivative set object
Definition xc_derivative_set_types.F:111
xc_derivatives
Definition xc_derivatives.F:9
xc_derivatives::xc_functionals_get_needs
type(xc_rho_cflags_type) function, public xc_functionals_get_needs(functionals, lsd, calc_potential)
...
Definition xc_derivatives.F:774
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::skala_gapw_direct_valence
integer, parameter, public skala_gapw_direct_valence
Definition xc_input_constants.F:29
xc_input_constants::skala_gapw_paw_one_center
integer, parameter, public skala_gapw_paw_one_center
Definition xc_input_constants.F:29
xc_input_constants::skala_gapw_cp2k_default
integer, parameter, public skala_gapw_cp2k_default
Definition xc_input_constants.F:29
xc_input_constants::skala_gapw_paw_one_center_split
integer, parameter, public skala_gapw_paw_one_center_split
Definition xc_input_constants.F:29
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
xc_rho_set_types
contains the structure
Definition xc_rho_set_types.F:14
xc_rho_set_types::xc_rho_set_create
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
Definition xc_rho_set_types.F:111
xc_rho_set_types::xc_rho_set_release
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
Definition xc_rho_set_types.F:129
xc_rho_set_types::xc_rho_set_update
subroutine, public xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, xc_deriv_method_id, xc_rho_smooth_id, pw_pool, spinflip)
updates the given rho set with the density given by rho_r (and rho_g). The rho set will contain the c...
Definition xc_rho_set_types.F:691
xc_rho_set_types::xc_rho_set_get
subroutine, public xc_rho_set_get(rho_set, can_return_null, rho, drho, norm_drho, rhoa, rhob, norm_drhoa, norm_drhob, rho_1_3, rhoa_1_3, rhob_1_3, laplace_rho, laplace_rhoa, laplace_rhob, drhoa, drhob, rho_cutoff, drho_cutoff, tau_cutoff, tau, tau_a, tau_b, local_bounds)
returns the various attributes of rho_set
Definition xc_rho_set_types.F:281
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_array_utils::cp_3d_r_cp_type
represent a pointer to a contiguous 3d array
Definition cp_array_utils.F:149
cp_control_types::dft_control_type
Definition cp_control_types.F:787
external_potential_types::gth_potential_type
Definition external_potential_types.F:121
external_potential_types::sgp_potential_type
Definition external_potential_types.F:179
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:53
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_cdft_grid::cdft_point_context_type
Definition qs_cdft_grid.F:41
qs_cdft_types::cdft_control_type
Definition qs_cdft_types.F:219
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_grid_atom::grid_atom_type
Definition qs_grid_atom.F:44
qs_harmonics_atom::harmonics_atom_type
Definition qs_harmonics_atom.F:32
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_linres_types::nablavks_atom_type
Definition qs_linres_types.F:213
qs_native_grid_cache::native_grid_interpolation_stencil_type
Definition qs_native_grid_cache.F:21
qs_rho_atom_types::rho_atom_coeff
Definition qs_rho_atom_types.F:19
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_vxc_atom_utils::tau_basis_cache_type
Definition qs_vxc_atom_utils.F:37
virial_types::virial_type
Definition virial_types.F:26
xc_derivative_set_types::xc_derivative_set_type
A derivative set contains the different derivatives of a xc-functional in form of a linked list.
Definition xc_derivative_set_types.F:48
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48
xc_rho_set_types::xc_rho_set_type
represent a density, with all the representation and data needed to perform a functional evaluation
Definition xc_rho_set_types.F:78