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