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