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