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 : #ifdef __GAUXC
9 : #include "gauxc/gauxc_config.f"
10 : #endif
11 :
12 : MODULE xc_gauxc_functional
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind
15 : USE cell_types, ONLY: cell_type
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
18 : dbcsr_create,&
19 : dbcsr_finalize,&
20 : dbcsr_get_info,&
21 : dbcsr_p_type,&
22 : dbcsr_release
23 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
24 : dbcsr_deallocate_matrix_set
25 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
26 : USE external_potential_types, ONLY: all_potential_type,&
27 : gth_potential_type,&
28 : sgp_potential_type
29 : USE input_constants, ONLY: xc_vdw_fun_nonloc
30 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
31 : section_vals_get_subs_vals2,&
32 : section_vals_type,&
33 : section_vals_val_get
34 : USE iso_c_binding, ONLY: c_char,&
35 : c_double,&
36 : c_int,&
37 : c_null_char
38 : USE kinds, ONLY: default_path_length,&
39 : default_string_length,&
40 : dp
41 : USE message_passing, ONLY: mp_comm_self,&
42 : mp_para_env_type
43 : USE particle_types, ONLY: particle_type
44 : USE qs_energy_types, ONLY: qs_energy_type
45 : USE qs_environment_types, ONLY: get_qs_env,&
46 : qs_environment_type
47 : USE qs_force_types, ONLY: qs_force_type
48 : USE qs_kind_types, ONLY: get_qs_kind,&
49 : has_nlcc,&
50 : qs_kind_type
51 : USE qs_ks_types, ONLY: qs_ks_env_type,&
52 : set_ks_env
53 : USE qs_rho_types, ONLY: qs_rho_get,&
54 : qs_rho_type
55 : USE qs_scf_types, ONLY: qs_scf_env_type
56 : USE string_utilities, ONLY: uppercase
57 : USE xc_gauxc_cache, ONLY: cp_gauxc_cache_params,&
58 : cp_gauxc_cache_type,&
59 : gauxc_cache_init
60 : USE xc_gauxc_interface, ONLY: &
61 : cp_gauxc_basisset_type, cp_gauxc_grid_type, cp_gauxc_integrator_type, &
62 : cp_gauxc_molecule_type, cp_gauxc_status_type, cp_gauxc_xc_gradient_type, cp_gauxc_xc_type, &
63 : gauxc_check_status, gauxc_compute_xc, gauxc_compute_xc_gradient, gauxc_create_basisset, &
64 : gauxc_create_grid, gauxc_create_integrator, gauxc_create_molecule, gauxc_destroy_basisset, &
65 : gauxc_destroy_grid, gauxc_destroy_integrator, gauxc_destroy_molecule, &
66 : gauxc_write_basisset_hdf5, gauxc_write_molecule_hdf5
67 : USE xc_input_constants, ONLY: skala_gapw_paw_one_center
68 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
69 : #include "../base/base_uses.f90"
70 :
71 : IMPLICIT NONE
72 :
73 : PRIVATE
74 :
75 : LOGICAL, PARAMETER :: debug_this_module = .TRUE.
76 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_gauxc_functional'
77 :
78 : PUBLIC :: apply_gauxc, gauxc_gapw_all_all_electron, gauxc_gapw_all_pseudopotentials, &
79 : gauxc_gapw_has_paw_pseudopotentials, gauxc_gapw_has_pseudopotentials, skala_info, &
80 : xc_section_uses_gauxc
81 :
82 : INTERFACE
83 : INTEGER(c_int) FUNCTION c_setenv(name, value, overwrite) BIND(C, name="setenv")
84 : IMPORT :: c_char, c_int
85 : CHARACTER(KIND=c_char), DIMENSION(*), INTENT(IN) :: name, value
86 : INTEGER(c_int), VALUE :: overwrite
87 : END FUNCTION c_setenv
88 :
89 : INTEGER(c_int) FUNCTION c_unsetenv(name) BIND(C, name="unsetenv")
90 : IMPORT :: c_char, c_int
91 : CHARACTER(KIND=c_char), DIMENSION(*), INTENT(IN) :: name
92 : END FUNCTION c_unsetenv
93 : END INTERFACE
94 :
95 : CONTAINS
96 :
97 : ! **************************************************************************************************
98 : !> \brief Return whether every GAPW kind uses an all-electron potential.
99 : !> \param qs_kind_set ...
100 : !> \return ...
101 : ! **************************************************************************************************
102 264 : FUNCTION gauxc_gapw_all_all_electron(qs_kind_set) RESULT(all_all_electron)
103 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
104 : LOGICAL :: all_all_electron
105 :
106 : INTEGER :: ikind
107 : TYPE(all_potential_type), POINTER :: all_potential
108 :
109 264 : CPASSERT(ASSOCIATED(qs_kind_set))
110 :
111 264 : all_all_electron = .TRUE.
112 384 : DO ikind = 1, SIZE(qs_kind_set)
113 272 : NULLIFY (all_potential)
114 272 : CALL get_qs_kind(qs_kind_set(ikind), all_potential=all_potential)
115 384 : IF (.NOT. ASSOCIATED(all_potential)) THEN
116 : all_all_electron = .FALSE.
117 : EXIT
118 : END IF
119 : END DO
120 :
121 264 : END FUNCTION gauxc_gapw_all_all_electron
122 :
123 : ! **************************************************************************************************
124 : !> \brief Set the GauXC Skala atom chunk environment knob when the CP2K keyword is explicit.
125 : !> \param atom_chunk_size ...
126 : !> \param is_explicit ...
127 : ! **************************************************************************************************
128 160 : SUBROUTINE set_gauxc_model_atom_chunk_env(atom_chunk_size, is_explicit)
129 : INTEGER, INTENT(IN) :: atom_chunk_size
130 : LOGICAL, INTENT(IN) :: is_explicit
131 :
132 : CHARACTER(LEN=32) :: chunk_value
133 : INTEGER(c_int) :: ierr
134 :
135 160 : IF (.NOT. is_explicit) RETURN
136 :
137 38 : IF (atom_chunk_size < 0) THEN
138 0 : ierr = c_unsetenv("GAUXC_ONEDFT_ATOM_CHUNK_SIZE"//c_null_char)
139 : ELSE
140 38 : WRITE (chunk_value, '(I0)') atom_chunk_size
141 : ierr = c_setenv( &
142 : "GAUXC_ONEDFT_ATOM_CHUNK_SIZE"//c_null_char, &
143 : TRIM(chunk_value)//c_null_char, &
144 38 : 1_c_int)
145 : END IF
146 38 : IF (ierr /= 0_c_int) THEN
147 : CALL cp_abort(__LOCATION__, &
148 0 : "Could not set GAUXC_ONEDFT_ATOM_CHUNK_SIZE for GauXC Skala.")
149 : END IF
150 : END SUBROUTINE set_gauxc_model_atom_chunk_env
151 :
152 : ! **************************************************************************************************
153 : !> \brief ...
154 : !> \param dbcsr_mat ...
155 : !> \param dense_mat ...
156 : !> \param para_env ...
157 : ! **************************************************************************************************
158 1576 : SUBROUTINE dbcsr_to_dense(dbcsr_mat, dense_mat, para_env)
159 : USE cp_dbcsr_api, ONLY: dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_info, &
160 : dbcsr_get_matrix_type, dbcsr_iterator_blocks_left, &
161 : dbcsr_iterator_next_block, dbcsr_iterator_readonly_start, &
162 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_type_antisymmetric, &
163 : dbcsr_type_symmetric
164 : TYPE(dbcsr_p_type), INTENT(IN) :: dbcsr_mat
165 : REAL(c_double), ALLOCATABLE, DIMENSION(:, :), &
166 : INTENT(INOUT) :: dense_mat
167 : TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
168 :
169 : CHARACTER :: matrix_type
170 : INTEGER :: col_end, col_start, icol, irow, ncol, &
171 : ncols, nrow, nrows, numnodes, row_end, &
172 : row_start
173 : LOGICAL :: has_symmetry, transposed
174 788 : REAL(c_double), POINTER :: block(:, :)
175 : REAL(KIND=dp) :: symmetry_sign
176 : TYPE(dbcsr_distribution_type) :: dist
177 : TYPE(dbcsr_iterator_type) :: iter
178 :
179 : CALL dbcsr_get_info(dbcsr_mat%matrix, &
180 : nfullrows_total=nrows, &
181 : nfullcols_total=ncols, &
182 788 : distribution=dist)
183 788 : CALL dbcsr_distribution_get(dist, numnodes=numnodes)
184 788 : matrix_type = dbcsr_get_matrix_type(dbcsr_mat%matrix)
185 788 : has_symmetry = matrix_type == dbcsr_type_symmetric .OR. matrix_type == dbcsr_type_antisymmetric
186 788 : symmetry_sign = 1.0_dp
187 788 : IF (matrix_type == dbcsr_type_antisymmetric) symmetry_sign = -1.0_dp
188 :
189 788 : IF (.NOT. ALLOCATED(dense_mat)) THEN
190 3152 : ALLOCATE (dense_mat(nrows, ncols))
191 0 : ELSE IF (.NOT. ALL(SHAPE(dense_mat) == [nrows, ncols])) THEN
192 0 : DEALLOCATE (dense_mat)
193 0 : ALLOCATE (dense_mat(nrows, ncols))
194 : ELSE
195 0 : CPASSERT(ALL(SHAPE(dense_mat) == [nrows, ncols]))
196 : END IF
197 788 : dense_mat = 0._dp
198 :
199 : ! Replicated DBCSR blocks must enter the following MPI sum exactly once.
200 788 : IF (numnodes == 1 .AND. para_env%num_pe > 1 .AND. para_env%mepos /= 0) RETURN
201 :
202 : ! Stored blocks and their symmetry mirrors occupy disjoint dense slices.
203 : !$OMP PARALLEL DEFAULT(NONE) SHARED(dbcsr_mat, dense_mat, has_symmetry, symmetry_sign) &
204 788 : !$OMP PRIVATE(iter, irow, icol, block, transposed, row_start, col_start, nrow, ncol, row_end, col_end)
205 : CALL dbcsr_iterator_readonly_start(iter, dbcsr_mat%matrix)
206 : DO WHILE (dbcsr_iterator_blocks_left(iter))
207 : CALL dbcsr_iterator_next_block(iter, irow, icol, block, transposed=transposed, &
208 : row_size=nrow, col_size=ncol, row_offset=row_start, col_offset=col_start)
209 : row_end = row_start + nrow - 1
210 : col_end = col_start + ncol - 1
211 : IF (transposed) THEN
212 : dense_mat(row_start:row_end, col_start:col_end) = symmetry_sign*TRANSPOSE(block)
213 : ELSE
214 : dense_mat(row_start:row_end, col_start:col_end) = block
215 : END IF
216 : IF (has_symmetry .AND. irow /= icol) THEN
217 : dense_mat(col_start:col_end, row_start:row_end) = &
218 : symmetry_sign*TRANSPOSE(dense_mat(row_start:row_end, col_start:col_end))
219 : END IF
220 : END DO
221 : CALL dbcsr_iterator_stop(iter)
222 : !$OMP END PARALLEL
223 :
224 788 : END SUBROUTINE dbcsr_to_dense
225 :
226 : ! ******, ***********************************************************************************
227 : !> \brief Convert a dense symmetric matrix to a DBCSR matrix with full upper block structure.
228 : !> This creates all upper-triangular blocks, not just those present in a template.
229 : !> This is needed because GauXC computes VXC for the full dense density matrix.
230 : !> \param dense_mat Input dense matrix
231 : !> \param template_dbcsr Template DBCSR matrix for distribution and block sizes
232 : !> \return dbcsr_mat Output DBCSR matrix with full upper block structure
233 : ! **************************************************************************************************
234 1812 : FUNCTION dense_to_dbcsr(dense_mat, template_dbcsr) RESULT(dbcsr_mat)
235 : USE cp_dbcsr_api, ONLY: &
236 : dbcsr_create, &
237 : dbcsr_distribution_get, &
238 : dbcsr_distribution_type, &
239 : dbcsr_finalize, &
240 : dbcsr_get_info, &
241 : dbcsr_get_stored_coordinates, &
242 : dbcsr_init_p, &
243 : dbcsr_put_block, &
244 : dbcsr_release, &
245 : dbcsr_type_symmetric, &
246 : dbcsr_work_create
247 : REAL(c_double), DIMENSION(:, :), INTENT(IN) :: dense_mat
248 : TYPE(dbcsr_p_type), INTENT(IN) :: template_dbcsr
249 : TYPE(dbcsr_p_type) :: dbcsr_mat
250 :
251 : INTEGER :: col, icol, irow, mynode, nblkcols_total, &
252 : nblkrows_total, ncols, nrows, owner, &
253 : row
254 906 : INTEGER, ALLOCATABLE, DIMENSION(:) :: c_offset, r_offset
255 906 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
256 : TYPE(dbcsr_distribution_type) :: dist
257 :
258 : CALL dbcsr_get_info(template_dbcsr%matrix, &
259 : row_blk_size=row_blk_size, &
260 : col_blk_size=col_blk_size, &
261 : nblkrows_total=nblkrows_total, &
262 : nblkcols_total=nblkcols_total, &
263 : nfullrows_total=nrows, &
264 : nfullcols_total=ncols, &
265 906 : distribution=dist)
266 906 : CALL dbcsr_distribution_get(dist, mynode=mynode)
267 :
268 906 : CPASSERT(nrows == SIZE(dense_mat, 1))
269 906 : CPASSERT(ncols == SIZE(dense_mat, 2))
270 :
271 906 : CALL dbcsr_init_p(dbcsr_mat%matrix)
272 : CALL dbcsr_create(dbcsr_mat%matrix, &
273 : template=template_dbcsr%matrix, &
274 : name="VXC from GauXC (dense)", &
275 906 : matrix_type=dbcsr_type_symmetric)
276 906 : CALL dbcsr_work_create(dbcsr_mat%matrix, work_mutable=.TRUE.)
277 :
278 4530 : ALLOCATE (r_offset(nblkrows_total), c_offset(nblkcols_total))
279 :
280 906 : r_offset(1) = 1
281 2028 : DO row = 2, nblkrows_total
282 2028 : r_offset(row) = r_offset(row - 1) + row_blk_size(row - 1)
283 : END DO
284 906 : c_offset(1) = 1
285 2028 : DO col = 2, nblkcols_total
286 2028 : c_offset(col) = c_offset(col - 1) + col_blk_size(col - 1)
287 : END DO
288 :
289 2934 : DO irow = 1, nblkrows_total
290 7950 : DO icol = 1, nblkcols_total
291 5016 : IF (irow > icol) CYCLE
292 3522 : CALL dbcsr_get_stored_coordinates(dbcsr_mat%matrix, irow, icol, owner)
293 3522 : IF (owner /= mynode) CYCLE
294 : CALL dbcsr_put_block(dbcsr_mat%matrix, irow, icol, &
295 : 0.5_dp*( &
296 : dense_mat(r_offset(irow):r_offset(irow) + row_blk_size(irow) - 1, &
297 : c_offset(icol):c_offset(icol) + col_blk_size(icol) - 1) + &
298 : TRANSPOSE(dense_mat(r_offset(icol):r_offset(icol) + row_blk_size(icol) - 1, &
299 165605 : c_offset(irow):c_offset(irow) + col_blk_size(irow) - 1))))
300 : END DO
301 : END DO
302 :
303 906 : CALL dbcsr_finalize(dbcsr_mat%matrix)
304 :
305 906 : DEALLOCATE (r_offset, c_offset)
306 :
307 906 : END FUNCTION dense_to_dbcsr
308 :
309 : ! **************************************************************************************************
310 : !> \brief ...
311 : !> \param xc_section ...
312 : !> \return ...
313 : ! **************************************************************************************************
314 670 : FUNCTION get_gauxc_functional(xc_section) RESULT(gauxc_functional_section)
315 : TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
316 : TYPE(section_vals_type), POINTER :: gauxc_functional_section
317 :
318 : INTEGER :: ifun
319 : TYPE(section_vals_type), POINTER :: functionals, xc_fun
320 :
321 670 : NULLIFY (gauxc_functional_section)
322 :
323 670 : functionals => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
324 670 : IF (.NOT. ASSOCIATED(functionals)) THEN
325 0 : CPABORT("XC_FUNCTIONAL section not found")
326 : END IF
327 :
328 670 : ifun = 0
329 : DO
330 1340 : ifun = ifun + 1
331 1340 : xc_fun => section_vals_get_subs_vals2(functionals, i_section=ifun)
332 1340 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
333 670 : IF (xc_fun%section%name /= "GAUXC" .OR. ifun > 1) THEN
334 0 : CPABORT("GauXC functionals are mutually exclusive with any other functional.")
335 : END IF
336 670 : gauxc_functional_section => xc_fun
337 : END DO
338 :
339 670 : IF (.NOT. ASSOCIATED(gauxc_functional_section)) THEN
340 0 : CPABORT("No XC functional found in XC_FUNCTIONAL section")
341 : END IF
342 670 : END FUNCTION get_gauxc_functional
343 :
344 : ! **************************************************************************************************
345 : !> \brief ...
346 : !> \param xc_section ...
347 : !> \return ...
348 : ! **************************************************************************************************
349 6902 : FUNCTION xc_section_uses_gauxc(xc_section) RESULT(uses_gauxc)
350 : TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
351 : LOGICAL :: uses_gauxc
352 :
353 : INTEGER :: ifun
354 : TYPE(section_vals_type), POINTER :: functionals, xc_fun
355 :
356 6902 : uses_gauxc = .FALSE.
357 6902 : IF (.NOT. ASSOCIATED(xc_section)) RETURN
358 :
359 6902 : functionals => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
360 6902 : IF (.NOT. ASSOCIATED(functionals)) RETURN
361 :
362 6902 : ifun = 0
363 : DO
364 13500 : ifun = ifun + 1
365 13500 : xc_fun => section_vals_get_subs_vals2(functionals, i_section=ifun)
366 13500 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
367 13500 : IF (xc_fun%section%name == "GAUXC") THEN
368 : uses_gauxc = .TRUE.
369 : EXIT
370 : END IF
371 : END DO
372 :
373 : END FUNCTION xc_section_uses_gauxc
374 :
375 : ! **************************************************************************************************
376 : !> \brief Return whether GauXC GAPW mode sees pseudopotential kinds.
377 : !> \param qs_kind_set ...
378 : !> \return ...
379 : ! **************************************************************************************************
380 206 : FUNCTION gauxc_gapw_has_pseudopotentials(qs_kind_set) RESULT(has_pseudopotentials)
381 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
382 : LOGICAL :: has_pseudopotentials
383 :
384 : INTEGER :: ikind
385 : TYPE(gth_potential_type), POINTER :: gth_potential
386 : TYPE(sgp_potential_type), POINTER :: sgp_potential
387 :
388 206 : CPASSERT(ASSOCIATED(qs_kind_set))
389 :
390 206 : has_pseudopotentials = .FALSE.
391 308 : DO ikind = 1, SIZE(qs_kind_set)
392 206 : NULLIFY (gth_potential, sgp_potential)
393 : CALL get_qs_kind(qs_kind_set(ikind), &
394 : gth_potential=gth_potential, &
395 206 : sgp_potential=sgp_potential)
396 308 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
397 : has_pseudopotentials = .TRUE.
398 : EXIT
399 : END IF
400 : END DO
401 :
402 206 : END FUNCTION gauxc_gapw_has_pseudopotentials
403 :
404 : ! **************************************************************************************************
405 : !> \brief Return whether every GAPW kind uses a pseudopotential.
406 : !> \param qs_kind_set ...
407 : !> \return ...
408 : ! **************************************************************************************************
409 168 : FUNCTION gauxc_gapw_all_pseudopotentials(qs_kind_set) RESULT(all_pseudopotentials)
410 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
411 : LOGICAL :: all_pseudopotentials
412 :
413 : INTEGER :: ikind
414 : TYPE(gth_potential_type), POINTER :: gth_potential
415 : TYPE(sgp_potential_type), POINTER :: sgp_potential
416 :
417 168 : CPASSERT(ASSOCIATED(qs_kind_set))
418 :
419 168 : all_pseudopotentials = .TRUE.
420 344 : DO ikind = 1, SIZE(qs_kind_set)
421 200 : NULLIFY (gth_potential, sgp_potential)
422 : CALL get_qs_kind(qs_kind_set(ikind), &
423 : gth_potential=gth_potential, &
424 200 : sgp_potential=sgp_potential)
425 344 : IF (.NOT. (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential))) THEN
426 : all_pseudopotentials = .FALSE.
427 : EXIT
428 : END IF
429 : END DO
430 :
431 168 : END FUNCTION gauxc_gapw_all_pseudopotentials
432 :
433 : ! **************************************************************************************************
434 : !> \brief Return whether GauXC GAPW mode sees pseudopotential one-center GAPW kinds.
435 : !> \param qs_kind_set ...
436 : !> \return ...
437 : ! **************************************************************************************************
438 492 : FUNCTION gauxc_gapw_has_paw_pseudopotentials(qs_kind_set) RESULT(has_paw_pseudopotentials)
439 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
440 : LOGICAL :: has_paw_pseudopotentials
441 :
442 : INTEGER :: ikind
443 : LOGICAL :: paw_atom
444 : TYPE(gth_potential_type), POINTER :: gth_potential
445 : TYPE(sgp_potential_type), POINTER :: sgp_potential
446 :
447 492 : CPASSERT(ASSOCIATED(qs_kind_set))
448 :
449 492 : has_paw_pseudopotentials = .FALSE.
450 1044 : DO ikind = 1, SIZE(qs_kind_set)
451 696 : NULLIFY (gth_potential, sgp_potential)
452 : CALL get_qs_kind(qs_kind_set(ikind), &
453 : gth_potential=gth_potential, &
454 : paw_atom=paw_atom, &
455 696 : sgp_potential=sgp_potential)
456 696 : IF ((ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) .AND. &
457 1044 : paw_atom) THEN
458 : has_paw_pseudopotentials = .TRUE.
459 : EXIT
460 : END IF
461 : END DO
462 :
463 492 : END FUNCTION gauxc_gapw_has_paw_pseudopotentials
464 :
465 : ! **************************************************************************************************
466 : !> \brief Check the current periodic scope of the CP2K-GauXC bridge
467 : !> \param dft_control ...
468 : !> \param cell ...
469 : !> \param qs_kind_set ...
470 : !> \param do_kpoints ...
471 : !> \param periodic_reference ...
472 : !> \note This path keeps isolated validation cells usable under PERIODIC XYZ.
473 : !> It intentionally does not implement compact periodic GauXC quadrature.
474 : ! **************************************************************************************************
475 670 : SUBROUTINE ensure_gauxc_periodic_reference_scope( &
476 : dft_control, cell, qs_kind_set, do_kpoints, periodic_reference)
477 : TYPE(dft_control_type), POINTER :: dft_control
478 : TYPE(cell_type), POINTER :: cell
479 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
480 : LOGICAL, INTENT(IN) :: do_kpoints, periodic_reference
481 :
482 : INTEGER :: ikind
483 : LOGICAL :: is_periodic
484 : TYPE(gth_potential_type), POINTER :: gth_potential
485 : TYPE(sgp_potential_type), POINTER :: sgp_potential
486 :
487 670 : CPASSERT(ASSOCIATED(dft_control))
488 670 : CPASSERT(ASSOCIATED(qs_kind_set))
489 :
490 670 : is_periodic = .FALSE.
491 1294 : IF (ASSOCIATED(cell)) is_periodic = ANY(cell%perd /= 0)
492 :
493 670 : IF (do_kpoints) THEN
494 : CALL cp_abort(__LOCATION__, &
495 : "GauXC currently supports only Gamma-only density matrices in CP2K. "// &
496 0 : "Periodic k-point density matrices require a dedicated GauXC periodic interface.")
497 : END IF
498 670 : IF (dft_control%nimages /= 1) THEN
499 : CALL cp_abort(__LOCATION__, &
500 : "GauXC currently supports only a single AO image in CP2K. "// &
501 0 : "Periodic neighbour-cell AO blocks require a dedicated GauXC periodic interface.")
502 : END IF
503 670 : IF (.NOT. is_periodic) RETURN
504 :
505 462 : IF (.NOT. periodic_reference) THEN
506 : CALL cp_abort(__LOCATION__, &
507 : "Periodic GauXC calculations in CP2K require GAUXC%PERIODIC_REFERENCE T. "// &
508 : "This opt-in documents that the current path is only an isolated-cell, "// &
509 : "Gamma-only, single-image METHOD GPW reference path using GauXC molecular "// &
510 0 : "quadrature, not a dedicated periodic GauXC interface.")
511 : END IF
512 :
513 1848 : IF (.NOT. ALL(cell%perd == 1)) THEN
514 : CALL cp_abort(__LOCATION__, &
515 : "The current GauXC isolated-cell reference path supports only PERIODIC XYZ. "// &
516 0 : "Partial periodicity requires a dedicated GauXC periodic interface.")
517 : END IF
518 462 : IF (.NOT. dft_control%qs_control%gpw) THEN
519 : CALL cp_abort(__LOCATION__, &
520 : "The current GauXC isolated-cell reference path is limited to METHOD GPW with GTH "// &
521 0 : "pseudopotentials. GAPW, GAPW_XC, and other QS methods are not supported here.")
522 : END IF
523 :
524 1038 : DO ikind = 1, SIZE(qs_kind_set)
525 576 : NULLIFY (gth_potential, sgp_potential)
526 : CALL get_qs_kind(qs_kind_set(ikind), &
527 : gth_potential=gth_potential, &
528 576 : sgp_potential=sgp_potential)
529 1038 : IF (.NOT. ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
530 : CALL cp_abort(__LOCATION__, &
531 : "The current GauXC isolated-cell reference path is limited to GTH pseudopotentials. "// &
532 0 : "Use non-periodic all-electron GAPW validation for molecular GAPW cases.")
533 : END IF
534 : END DO
535 :
536 : END SUBROUTINE ensure_gauxc_periodic_reference_scope
537 :
538 : ! **************************************************************************************************
539 : !> \brief adds a replicated GauXC energy gradient to the local CP2K force accumulator
540 : !> \param exc_grad ...
541 : !> \param force ...
542 : !> \param atomic_kind_set ...
543 : !> \param para_env ...
544 : ! **************************************************************************************************
545 10 : SUBROUTINE add_gauxc_gradient_to_force(exc_grad, force, atomic_kind_set, para_env)
546 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
547 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
548 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
549 : TYPE(mp_para_env_type), POINTER :: para_env
550 :
551 : INTEGER :: ia, iatom, ikind, natom_kind
552 : TYPE(atomic_kind_type), POINTER :: atomic_kind
553 :
554 10 : CPASSERT(ASSOCIATED(force))
555 10 : CPASSERT(ASSOCIATED(atomic_kind_set))
556 :
557 10 : IF (para_env%mepos /= 0) RETURN
558 :
559 12 : DO ikind = 1, SIZE(atomic_kind_set, 1)
560 7 : atomic_kind => atomic_kind_set(ikind)
561 7 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
562 24 : DO ia = 1, natom_kind
563 12 : iatom = atomic_kind%atom_list(ia)
564 : force(ikind)%rho_elec(:, ia) = force(ikind)%rho_elec(:, ia) + &
565 55 : exc_grad(3*iatom - 2:3*iatom)
566 : END DO
567 : END DO
568 :
569 : END SUBROUTINE add_gauxc_gradient_to_force
570 :
571 : ! **************************************************************************************************
572 : !> \brief compute a GauXC XC energy for diagnostic finite differences
573 : !> \param particle_set_eval ...
574 : !> \param qs_kind_set ...
575 : !> \param density_scalar ...
576 : !> \param nspins ...
577 : !> \param model_name ...
578 : !> \param xc_fun_name ...
579 : !> \param grid_type ...
580 : !> \param radial_quadrature ...
581 : !> \param pruning_scheme ...
582 : !> \param lb_exec_space ...
583 : !> \param int_exec_space ...
584 : !> \param lwd_kernel ...
585 : !> \param batch_size ...
586 : !> \param device_runtime_fill_fraction ...
587 : !> \param exc ...
588 : !> \param density_zeta ...
589 : ! **************************************************************************************************
590 60 : SUBROUTINE gauxc_xc_energy_for_particles( &
591 60 : particle_set_eval, qs_kind_set, density_scalar, nspins, model_name, &
592 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
593 60 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, exc, density_zeta)
594 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_eval
595 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
596 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
597 : INTEGER, INTENT(IN) :: nspins
598 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
599 : pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
600 : INTEGER, INTENT(IN) :: batch_size
601 : REAL(KIND=dp), INTENT(IN) :: device_runtime_fill_fraction
602 : REAL(KIND=dp), INTENT(OUT) :: exc
603 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
604 : OPTIONAL :: density_zeta
605 :
606 : TYPE(cp_gauxc_basisset_type) :: gauxc_basis_fd
607 : TYPE(cp_gauxc_grid_type) :: gauxc_grid_fd
608 : TYPE(cp_gauxc_integrator_type) :: gauxc_integrator_fd
609 : TYPE(cp_gauxc_molecule_type) :: gauxc_mol_fd
610 : TYPE(cp_gauxc_status_type) :: gauxc_status
611 60 : TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
612 :
613 60 : gauxc_mol_fd = gauxc_create_molecule(particle_set_eval, gauxc_status)
614 60 : CALL gauxc_check_status(gauxc_status)
615 60 : gauxc_basis_fd = gauxc_create_basisset(qs_kind_set, particle_set_eval, gauxc_status)
616 60 : CALL gauxc_check_status(gauxc_status)
617 : gauxc_grid_fd = gauxc_create_grid( &
618 : gauxc_mol_fd, &
619 : gauxc_basis_fd, &
620 : grid_type, &
621 : radial_quadrature, &
622 : pruning_scheme, &
623 : lb_exec_space, &
624 : batch_size, &
625 : device_runtime_fill_fraction, &
626 : gauxc_status, &
627 : mpi_comm=mp_comm_self%get_handle(), &
628 60 : force_new_runtime=.TRUE.)
629 60 : CALL gauxc_check_status(gauxc_status)
630 : gauxc_integrator_fd = gauxc_create_integrator( &
631 : TRIM(xc_fun_name), &
632 : gauxc_grid_fd, &
633 : int_exec_space, &
634 : lwd_kernel, &
635 : nspins, &
636 60 : gauxc_status)
637 60 : CALL gauxc_check_status(gauxc_status)
638 :
639 60 : IF (nspins == 1) THEN
640 : gauxc_xc_result = gauxc_compute_xc( &
641 : gauxc_integrator_fd, &
642 : density_scalar, &
643 : nspins=nspins, &
644 : status=gauxc_status, &
645 60 : model=TRIM(model_name))
646 : ELSE
647 0 : CPASSERT(nspins == 2)
648 0 : CPASSERT(PRESENT(density_zeta))
649 : gauxc_xc_result = gauxc_compute_xc( &
650 : gauxc_integrator_fd, &
651 : density_scalar, &
652 : density_zeta, &
653 : nspins, &
654 : gauxc_status, &
655 0 : model=TRIM(model_name))
656 : END IF
657 60 : CALL gauxc_check_status(gauxc_status)
658 60 : exc = gauxc_xc_result%exc
659 :
660 60 : IF (ALLOCATED(gauxc_xc_result%vxc_scalar)) DEALLOCATE (gauxc_xc_result%vxc_scalar)
661 60 : IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
662 :
663 60 : CALL gauxc_destroy_integrator(gauxc_integrator_fd, gauxc_status)
664 60 : CALL gauxc_check_status(gauxc_status)
665 60 : CALL gauxc_destroy_grid(gauxc_grid_fd, gauxc_status)
666 60 : CALL gauxc_check_status(gauxc_status)
667 60 : CALL gauxc_destroy_basisset(gauxc_basis_fd, gauxc_status)
668 60 : CALL gauxc_check_status(gauxc_status)
669 60 : CALL gauxc_destroy_molecule(gauxc_mol_fd, gauxc_status)
670 60 : CALL gauxc_check_status(gauxc_status)
671 :
672 120 : END SUBROUTINE gauxc_xc_energy_for_particles
673 :
674 : ! **************************************************************************************************
675 : !> \brief compute a finite-difference GauXC XC nuclear gradient at fixed density
676 : !> \param particle_set ...
677 : !> \param qs_kind_set ...
678 : !> \param density_scalar ...
679 : !> \param nspins ...
680 : !> \param model_name ...
681 : !> \param xc_fun_name ...
682 : !> \param grid_type ...
683 : !> \param radial_quadrature ...
684 : !> \param pruning_scheme ...
685 : !> \param lb_exec_space ...
686 : !> \param int_exec_space ...
687 : !> \param lwd_kernel ...
688 : !> \param batch_size ...
689 : !> \param device_runtime_fill_fraction ...
690 : !> \param dx ...
691 : !> \param para_env ...
692 : !> \param exc_grad ...
693 : !> \param density_zeta ...
694 : ! **************************************************************************************************
695 2 : SUBROUTINE gauxc_xc_gradient_fd( &
696 2 : particle_set, qs_kind_set, density_scalar, nspins, model_name, &
697 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
698 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, dx, para_env, exc_grad, &
699 2 : density_zeta)
700 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
701 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
702 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
703 : INTEGER, INTENT(IN) :: nspins
704 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
705 : pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
706 : INTEGER, INTENT(IN) :: batch_size
707 : REAL(KIND=dp), INTENT(IN) :: device_runtime_fill_fraction, dx
708 : TYPE(mp_para_env_type), POINTER :: para_env
709 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
710 : INTENT(OUT) :: exc_grad
711 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
712 : OPTIONAL :: density_zeta
713 :
714 : INTEGER :: iatom, idir
715 : REAL(KIND=dp) :: xc_minus, xc_plus
716 2 : TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus
717 :
718 2 : CPASSERT(ASSOCIATED(particle_set))
719 2 : CPASSERT(dx > 0.0_dp)
720 :
721 6 : ALLOCATE (exc_grad(3*SIZE(particle_set)))
722 2 : exc_grad = 0.0_dp
723 :
724 2 : IF (para_env%mepos == 0) THEN
725 31 : ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set)))
726 :
727 3 : DO iatom = 1, SIZE(particle_set)
728 9 : DO idir = 1, 3
729 24 : particle_set_minus = particle_set
730 24 : particle_set_plus = particle_set
731 6 : particle_set_minus(iatom)%r(idir) = particle_set_minus(iatom)%r(idir) - dx
732 6 : particle_set_plus(iatom)%r(idir) = particle_set_plus(iatom)%r(idir) + dx
733 6 : IF (PRESENT(density_zeta)) THEN
734 : CALL gauxc_xc_energy_for_particles( &
735 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
736 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
737 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus, &
738 0 : density_zeta=density_zeta)
739 : CALL gauxc_xc_energy_for_particles( &
740 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
741 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
742 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus, &
743 0 : density_zeta=density_zeta)
744 : ELSE
745 : CALL gauxc_xc_energy_for_particles( &
746 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
747 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
748 6 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus)
749 : CALL gauxc_xc_energy_for_particles( &
750 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
751 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
752 6 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus)
753 : END IF
754 8 : exc_grad(3*iatom - 3 + idir) = (xc_plus - xc_minus)/(2.0_dp*dx)
755 : END DO
756 : END DO
757 :
758 1 : DEALLOCATE (particle_set_minus, particle_set_plus)
759 : END IF
760 :
761 2 : CALL para_env%bcast(exc_grad, 0)
762 :
763 2 : END SUBROUTINE gauxc_xc_gradient_fd
764 :
765 : ! **************************************************************************************************
766 : !> \brief finite-difference check of the molecular GauXC XC virial diagnostic
767 : !> \param exc_grad ...
768 : !> \param particle_set ...
769 : !> \param qs_kind_set ...
770 : !> \param density_scalar ...
771 : !> \param nspins ...
772 : !> \param model_name ...
773 : !> \param xc_fun_name ...
774 : !> \param grid_type ...
775 : !> \param radial_quadrature ...
776 : !> \param pruning_scheme ...
777 : !> \param lb_exec_space ...
778 : !> \param int_exec_space ...
779 : !> \param lwd_kernel ...
780 : !> \param batch_size ...
781 : !> \param device_runtime_fill_fraction ...
782 : !> \param dx ...
783 : !> \param para_env ...
784 : !> \param density_zeta ...
785 : ! **************************************************************************************************
786 48 : SUBROUTINE debug_gauxc_molecular_virial( &
787 96 : exc_grad, particle_set, qs_kind_set, density_scalar, nspins, model_name, &
788 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
789 48 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, dx, para_env, density_zeta)
790 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
791 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
792 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
793 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
794 : INTEGER, INTENT(IN) :: nspins
795 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
796 : pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
797 : INTEGER, INTENT(IN) :: batch_size
798 : REAL(KIND=dp), INTENT(IN) :: device_runtime_fill_fraction, dx
799 : TYPE(mp_para_env_type), POINTER :: para_env
800 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
801 : OPTIONAL :: density_zeta
802 :
803 : INTEGER :: iatom, iw
804 : REAL(KIND=dp) :: analytic_trace, diff_trace, &
805 : numerical_trace, xc_minus, xc_plus
806 : REAL(KIND=dp), DIMENSION(3) :: center, displacement, grad
807 48 : TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus
808 :
809 48 : CPASSERT(ASSOCIATED(particle_set))
810 48 : CPASSERT(SIZE(exc_grad) == 3*SIZE(particle_set))
811 :
812 48 : IF (para_env%mepos /= 0) RETURN
813 :
814 24 : center = 0.0_dp
815 72 : DO iatom = 1, SIZE(particle_set)
816 216 : center = center + particle_set(iatom)%r
817 : END DO
818 96 : center = center/REAL(SIZE(particle_set), dp)
819 :
820 648 : ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set)))
821 96 : particle_set_minus = particle_set
822 96 : particle_set_plus = particle_set
823 :
824 24 : analytic_trace = 0.0_dp
825 72 : DO iatom = 1, SIZE(particle_set)
826 192 : grad = exc_grad(3*iatom - 2:3*iatom)
827 192 : displacement = particle_set(iatom)%r - center
828 192 : analytic_trace = analytic_trace + DOT_PRODUCT(grad, displacement)
829 192 : particle_set_minus(iatom)%r = center + (1.0_dp - dx)*displacement
830 216 : particle_set_plus(iatom)%r = center + (1.0_dp + dx)*displacement
831 : END DO
832 24 : analytic_trace = analytic_trace/3.0_dp
833 :
834 24 : IF (PRESENT(density_zeta)) THEN
835 : CALL gauxc_xc_energy_for_particles( &
836 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
837 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
838 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus, &
839 0 : density_zeta=density_zeta)
840 : CALL gauxc_xc_energy_for_particles( &
841 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
842 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
843 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus, &
844 0 : density_zeta=density_zeta)
845 : ELSE
846 : CALL gauxc_xc_energy_for_particles( &
847 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
848 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
849 24 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus)
850 : CALL gauxc_xc_energy_for_particles( &
851 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
852 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
853 24 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus)
854 : END IF
855 :
856 24 : numerical_trace = (xc_plus - xc_minus)/(2.0_dp*dx)/3.0_dp
857 24 : diff_trace = analytic_trace - numerical_trace
858 :
859 24 : iw = cp_logger_get_default_io_unit()
860 24 : IF (iw > 0) THEN
861 : WRITE (UNIT=iw, FMT="(/,T2,A,1X,ES11.4)") &
862 24 : "GAUXC| Molecular XC virial finite-difference dx", dx
863 : WRITE (UNIT=iw, FMT="(T2,A,3(1X,ES19.11))") &
864 24 : "GAUXC| Molecular XC virial FD 1/3 Trace", &
865 48 : analytic_trace, numerical_trace, diff_trace
866 : END IF
867 :
868 24 : DEALLOCATE (particle_set_minus, particle_set_plus)
869 :
870 : END SUBROUTINE debug_gauxc_molecular_virial
871 :
872 : ! **************************************************************************************************
873 : !> \brief prints a force-based molecular XC virial diagnostic from GauXC gradients
874 : !> \param exc_grad ...
875 : !> \param particle_set ...
876 : !> \param para_env ...
877 : ! **************************************************************************************************
878 48 : SUBROUTINE print_gauxc_molecular_virial(exc_grad, particle_set, para_env)
879 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
880 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
881 : TYPE(mp_para_env_type), POINTER :: para_env
882 :
883 : CHARACTER(len=1), DIMENSION(3), PARAMETER :: label = ["x", "y", "z"]
884 :
885 : INTEGER :: i, iatom, iw, j
886 : REAL(KIND=dp), DIMENSION(3) :: center, displacement, grad, grad_sum
887 : REAL(KIND=dp), DIMENSION(3, 3) :: molecular_virial
888 :
889 48 : CPASSERT(ASSOCIATED(particle_set))
890 48 : CPASSERT(SIZE(exc_grad) == 3*SIZE(particle_set))
891 :
892 48 : IF (para_env%mepos /= 0) RETURN
893 :
894 24 : center = 0.0_dp
895 72 : DO iatom = 1, SIZE(particle_set)
896 216 : center = center + particle_set(iatom)%r
897 : END DO
898 96 : center = center/REAL(SIZE(particle_set), dp)
899 :
900 24 : grad_sum = 0.0_dp
901 24 : molecular_virial = 0.0_dp
902 72 : DO iatom = 1, SIZE(particle_set)
903 192 : grad = exc_grad(3*iatom - 2:3*iatom)
904 192 : displacement = particle_set(iatom)%r - center
905 192 : grad_sum = grad_sum + grad
906 216 : DO i = 1, 3
907 624 : DO j = 1, 3
908 576 : molecular_virial(i, j) = molecular_virial(i, j) + grad(i)*displacement(j)
909 : END DO
910 : END DO
911 : END DO
912 :
913 24 : iw = cp_logger_get_default_io_unit()
914 24 : IF (iw <= 0) RETURN
915 :
916 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
917 24 : "GAUXC| Molecular XC gradient virial diagnostic [a.u.]"
918 24 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,T40,A,T60,A)") "GAUXC|", "x", "y", "z"
919 96 : DO i = 1, 3
920 : WRITE (UNIT=iw, FMT="(T2,A,1X,A1,3(1X,ES19.11))") &
921 96 : "GAUXC|", label(i), molecular_virial(i, :)
922 : END DO
923 : WRITE (UNIT=iw, FMT="(T2,A,1X,ES19.11)") &
924 24 : "GAUXC| Molecular XC gradient virial 1/3 Trace", &
925 48 : (molecular_virial(1, 1) + molecular_virial(2, 2) + molecular_virial(3, 3))/3.0_dp
926 : WRITE (UNIT=iw, FMT="(T2,A,3(1X,ES19.11))") &
927 24 : "GAUXC| Molecular XC gradient sum", grad_sum
928 : WRITE (UNIT=iw, FMT="(T2,A)") &
929 24 : "GAUXC| Diagnostic only; this is not an analytical periodic stress tensor."
930 :
931 : END SUBROUTINE print_gauxc_molecular_virial
932 :
933 : ! **************************************************************************************************
934 : !> \brief Return information about the Skala functional
935 : !> \param functional section containing the SKALA subsection
936 : !> \param lsd if you are using lsd or lda
937 : !> \param reference the reference to the article where the functional is explained
938 : !> \param shortform the short definition of the functional
939 : !> \param needs the flags corresponding to the inputs needed by this
940 : !> functional are set to true (the flags not needed aren't touched)
941 : !> \param max_deriv the maximal derivative available
942 : ! **************************************************************************************************
943 1557 : SUBROUTINE skala_info(functional, lsd, reference, shortform, needs, max_deriv)
944 : TYPE(section_vals_type), POINTER :: functional
945 : LOGICAL, INTENT(in) :: lsd
946 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
947 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
948 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
949 :
950 : CHARACTER(len=default_path_length) :: model_key, model_name
951 : CHARACTER(len=default_string_length) :: xc_fun_key, xc_fun_name
952 : INTEGER :: gapw_representation
953 : LOGICAL :: native_grid
954 :
955 519 : CALL section_vals_val_get(functional, "FUNCTIONAL", c_val=xc_fun_name)
956 519 : CALL section_vals_val_get(functional, "MODEL", c_val=model_name)
957 519 : CALL section_vals_val_get(functional, "NATIVE_GRID", l_val=native_grid)
958 : CALL section_vals_val_get(functional, "PSEUDOPOTENTIAL_GAPW_REPRESENTATION", &
959 519 : i_val=gapw_representation)
960 519 : native_grid = native_grid .OR. gapw_representation == skala_gapw_paw_one_center
961 519 : model_key = ADJUSTL(model_name)
962 519 : xc_fun_key = ADJUSTL(xc_fun_name)
963 519 : CALL uppercase(model_key)
964 519 : CALL uppercase(xc_fun_key)
965 :
966 519 : IF (PRESENT(reference)) THEN
967 5 : IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "" .OR. &
968 : TRIM(model_key) == TRIM(xc_fun_key)) THEN
969 0 : reference = "Functional computed by GauXC (underlying: "//TRIM(xc_fun_name)//")"
970 : ELSE
971 5 : reference = "Functional computed by GauXC Skala model "//TRIM(model_name)
972 : END IF
973 : END IF
974 519 : IF (PRESENT(shortform)) THEN
975 5 : IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "" .OR. &
976 : TRIM(model_key) == TRIM(xc_fun_key)) THEN
977 0 : shortform = "GAUXC ("//TRIM(xc_fun_name)//")"
978 : ELSE
979 5 : shortform = "GAUXC Skala"
980 : END IF
981 : END IF
982 519 : IF (PRESENT(needs)) THEN
983 514 : IF (native_grid .AND. TRIM(model_key) /= "NONE" .AND. TRIM(model_key) /= "" .AND. &
984 : TRIM(model_key) /= TRIM(xc_fun_key)) THEN
985 460 : IF (lsd) THEN
986 58 : needs%rho_spin = .TRUE.
987 58 : needs%drho_spin = .TRUE.
988 58 : needs%tau_spin = .TRUE.
989 : ELSE
990 402 : needs%rho = .TRUE.
991 402 : needs%drho = .TRUE.
992 402 : needs%tau = .TRUE.
993 : END IF
994 : ELSE
995 54 : needs%rho = .TRUE.
996 54 : IF (lsd) THEN
997 16 : needs%rho_spin = .TRUE.
998 : END IF
999 : END IF
1000 : END IF
1001 519 : IF (PRESENT(max_deriv)) max_deriv = 1
1002 :
1003 519 : END SUBROUTINE skala_info
1004 :
1005 : ! GauXC uses replicated dense density and VXC matrices. The DBCSR density matrix
1006 : ! is distributed over MPI ranks, so apply_gauxc allreduces the dense copy before
1007 : ! passing it to GauXC.
1008 :
1009 : ! **************************************************************************************************
1010 : !> \brief Evaluate molecular GauXC, optionally for an explicit configuration density.
1011 : !> \param qs_env ...
1012 : !> \param xc_section ...
1013 : !> \param calculate_forces ...
1014 : !> \param rho_ao_external Optional spin/image AO density; supply all three external arguments together
1015 : !> \param exc_external Configuration XC energy, leaving the baseline energy unchanged
1016 : !> \param matrix_vxc_external Initially null; caller owns the returned spin AO potential matrices
1017 : !> \note External-density mode preserves baseline energy and KS storage and shares
1018 : !> the geometry/model cache.
1019 : ! **************************************************************************************************
1020 670 : SUBROUTINE apply_gauxc(qs_env, xc_section, calculate_forces, &
1021 : rho_ao_external, exc_external, matrix_vxc_external)
1022 : TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
1023 : TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
1024 : LOGICAL, INTENT(IN) :: calculate_forces
1025 : TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
1026 : OPTIONAL, POINTER :: rho_ao_external
1027 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: exc_external
1028 : TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
1029 : OPTIONAL, POINTER :: matrix_vxc_external
1030 :
1031 : CHARACTER(len=*), PARAMETER :: nonlocal_vdw_abort_message = &
1032 : "GauXC does not support non-local VDW_POTENTIAL corrections. "// &
1033 : "Use an additive PAIR_POTENTIAL dispersion correction or disable GauXC."
1034 : REAL(KIND=dp), PARAMETER :: gapw_fd_gradient_dx = 1.0E-4_dp
1035 :
1036 : CHARACTER(len=default_path_length) :: model_key, model_name, output_path
1037 : CHARACTER(len=default_string_length) :: gradient_runtime, gradient_runtime_key, grid_key, &
1038 : int_exec_space_key, lwd_kernel_key, pruning_key, skala_runtime, skala_runtime_key, &
1039 : xc_fun_key
1040 : INTEGER :: atom_chunk_size, env_status, img, ispin, &
1041 : nimages
1042 : LOGICAL :: atom_chunk_size_explicit, do_kpoints, external_density, gapw_method, &
1043 : gapw_paw_pseudopotentials, gapw_pseudopotentials, grid_explicit, hdf5_output, &
1044 : is_periodic, low_spin_scope, molecular_virial, molecular_virial_debug, need_xc_gradient, &
1045 : periodic_reference, pruning_explicit, use_skala_model, write_hdf5_output
1046 : REAL(KIND=dp) :: exc_value, molecular_virial_debug_dx
1047 670 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: density_scalar, density_zeta
1048 670 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1049 : TYPE(cell_type), POINTER :: cell
1050 : TYPE(cp_gauxc_cache_params) :: params
1051 : TYPE(cp_gauxc_cache_type), POINTER :: cache
1052 : TYPE(cp_gauxc_status_type) :: gauxc_status
1053 670 : TYPE(cp_gauxc_xc_gradient_type) :: exc_grad
1054 670 : TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
1055 : TYPE(dbcsr_p_type) :: vxc_zeta_tmp
1056 670 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
1057 670 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
1058 : TYPE(dft_control_type), POINTER :: dft_control
1059 : TYPE(mp_para_env_type), POINTER :: para_env
1060 670 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1061 : TYPE(qs_energy_type), POINTER :: energy
1062 670 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1063 670 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1064 : TYPE(qs_ks_env_type), POINTER :: ks_env
1065 : TYPE(qs_rho_type), POINTER :: rho, rho_use, rho_xc
1066 : TYPE(qs_scf_env_type), POINTER :: scf_env
1067 : TYPE(section_vals_type), POINTER :: gauxc_functional_section
1068 :
1069 670 : external_density = PRESENT(rho_ao_external)
1070 0 : CPASSERT(PRESENT(exc_external) .EQV. external_density)
1071 670 : CPASSERT(PRESENT(matrix_vxc_external) .EQV. external_density)
1072 670 : IF (external_density) THEN
1073 64 : CPASSERT(ASSOCIATED(rho_ao_external))
1074 64 : CPASSERT(.NOT. ASSOCIATED(matrix_vxc_external))
1075 : END IF
1076 :
1077 : NULLIFY ( &
1078 : atomic_kind_set, &
1079 670 : cell, &
1080 670 : dft_control, &
1081 670 : energy, &
1082 670 : force, &
1083 670 : ks_env, &
1084 670 : matrix_vxc, &
1085 670 : para_env, &
1086 670 : particle_set, &
1087 670 : qs_kind_set, &
1088 670 : rho, &
1089 670 : rho_use, &
1090 670 : rho_xc, &
1091 670 : rho_ao, &
1092 670 : scf_env)
1093 :
1094 : CALL get_qs_env( &
1095 : qs_env, &
1096 : cell=cell, &
1097 : dft_control=dft_control, &
1098 : do_kpoints=do_kpoints, &
1099 : energy=energy, &
1100 : ks_env=ks_env, &
1101 : matrix_vxc=matrix_vxc, &
1102 : natom=params%natom, &
1103 : atomic_kind_set=atomic_kind_set, &
1104 : force=force, &
1105 : para_env=para_env, &
1106 : particle_set=particle_set, &
1107 : qs_kind_set=qs_kind_set, &
1108 : rho=rho, &
1109 : rho_xc=rho_xc, &
1110 670 : scf_env=scf_env)
1111 :
1112 670 : gapw_method = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1113 : gapw_pseudopotentials = gapw_method .AND. &
1114 206 : gauxc_gapw_has_pseudopotentials(qs_kind_set)
1115 : gapw_paw_pseudopotentials = gapw_method .AND. &
1116 206 : gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
1117 670 : IF (dft_control%qs_control%gapw_xc) THEN
1118 0 : CPASSERT(ASSOCIATED(rho_xc))
1119 0 : rho_use => rho_xc
1120 : ELSE
1121 670 : CPASSERT(ASSOCIATED(rho))
1122 670 : rho_use => rho
1123 : END IF
1124 : CALL qs_rho_get( &
1125 : rho_use, &
1126 670 : rho_ao_kp=rho_ao)
1127 :
1128 670 : IF (external_density) THEN
1129 64 : rho_ao => rho_ao_external
1130 : ! get_qs_env returned the baseline pointer. Never release or replace it.
1131 64 : NULLIFY (matrix_vxc)
1132 : END IF
1133 670 : nimages = dft_control%nimages
1134 670 : params%nspins = dft_control%nspins
1135 670 : IF (external_density) THEN
1136 64 : CPASSERT(SIZE(rho_ao, 1) == params%nspins)
1137 64 : CPASSERT(SIZE(rho_ao, 2) == nimages)
1138 : END IF
1139 606 : low_spin_scope = external_density .OR. (dft_control%low_spin_roks .AND. gapw_method)
1140 670 : is_periodic = .FALSE.
1141 1294 : IF (ASSOCIATED(cell)) is_periodic = ANY(cell%perd /= 0)
1142 :
1143 670 : IF (ASSOCIATED(qs_env%dispersion_env)) THEN
1144 670 : IF (qs_env%dispersion_env%type == xc_vdw_fun_nonloc) THEN
1145 0 : CPABORT(nonlocal_vdw_abort_message)
1146 : END IF
1147 : END IF
1148 670 : NULLIFY (vxc_zeta_tmp%matrix)
1149 :
1150 670 : gauxc_functional_section => get_gauxc_functional(xc_section)
1151 : CALL section_vals_val_get( &
1152 : gauxc_functional_section, &
1153 : "FUNCTIONAL", &
1154 670 : c_val=params%xc_fun_name)
1155 : CALL section_vals_val_get( &
1156 : gauxc_functional_section, &
1157 : "MODEL", &
1158 670 : c_val=model_name)
1159 : CALL section_vals_val_get( &
1160 : gauxc_functional_section, &
1161 : "GRID", &
1162 : c_val=params%grid_type, &
1163 670 : explicit=grid_explicit)
1164 : CALL section_vals_val_get( &
1165 : gauxc_functional_section, &
1166 : "RADIAL_QUADRATURE", &
1167 670 : c_val=params%radial_quadrature)
1168 : CALL section_vals_val_get( &
1169 : gauxc_functional_section, &
1170 : "PRUNING_SCHEME", &
1171 : c_val=params%pruning_scheme, &
1172 670 : explicit=pruning_explicit)
1173 : CALL section_vals_val_get( &
1174 : gauxc_functional_section, &
1175 : "BATCH_SIZE", &
1176 670 : i_val=params%batch_size)
1177 : CALL section_vals_val_get( &
1178 : gauxc_functional_section, &
1179 : "DEVICE_RUNTIME_FILL_FRACTION", &
1180 670 : r_val=params%device_runtime_fill_fraction)
1181 : CALL section_vals_val_get( &
1182 : gauxc_functional_section, &
1183 : "MODEL_ATOM_CHUNK_SIZE", &
1184 : i_val=atom_chunk_size, &
1185 670 : explicit=atom_chunk_size_explicit)
1186 : CALL section_vals_val_get( &
1187 : gauxc_functional_section, &
1188 : "PERIODIC_REFERENCE", &
1189 670 : l_val=periodic_reference)
1190 : CALL section_vals_val_get( &
1191 : gauxc_functional_section, &
1192 : "MOLECULAR_VIRIAL", &
1193 670 : l_val=molecular_virial)
1194 : CALL section_vals_val_get( &
1195 : gauxc_functional_section, &
1196 : "MOLECULAR_VIRIAL_DEBUG", &
1197 670 : l_val=molecular_virial_debug)
1198 : CALL section_vals_val_get( &
1199 : gauxc_functional_section, &
1200 : "MOLECULAR_VIRIAL_DEBUG_DX", &
1201 670 : r_val=molecular_virial_debug_dx)
1202 : CALL section_vals_val_get( &
1203 : gauxc_functional_section, &
1204 : "LB_EXECUTION_SPACE", &
1205 670 : c_val=params%lb_exec_space)
1206 : CALL section_vals_val_get( &
1207 : gauxc_functional_section, &
1208 : "INT_EXECUTION_SPACE", &
1209 670 : c_val=params%int_exec_space)
1210 : CALL section_vals_val_get( &
1211 : gauxc_functional_section, &
1212 : "LWD_KERNEL", &
1213 670 : c_val=params%lwd_kernel)
1214 : CALL section_vals_val_get( &
1215 : gauxc_functional_section, &
1216 : "SKALA_RUNTIME", &
1217 670 : c_val=skala_runtime)
1218 : CALL section_vals_val_get( &
1219 : gauxc_functional_section, &
1220 : "MODEL_GRADIENT_RUNTIME", &
1221 670 : c_val=gradient_runtime)
1222 : CALL section_vals_val_get( &
1223 : gauxc_functional_section, &
1224 : "FORCE_NUMERICAL_GRADIENT", &
1225 670 : l_val=params%force_fd_gradient)
1226 : CALL section_vals_val_get( &
1227 : gauxc_functional_section, &
1228 : "OUTPUT_PATH", &
1229 670 : c_val=output_path)
1230 :
1231 670 : model_key = ADJUSTL(model_name)
1232 670 : CALL uppercase(model_key)
1233 670 : xc_fun_key = ADJUSTL(params%xc_fun_name)
1234 670 : CALL uppercase(xc_fun_key)
1235 670 : skala_runtime_key = ADJUSTL(skala_runtime)
1236 670 : CALL uppercase(skala_runtime_key)
1237 670 : gradient_runtime_key = ADJUSTL(gradient_runtime)
1238 670 : CALL uppercase(gradient_runtime_key)
1239 670 : int_exec_space_key = ADJUSTL(params%int_exec_space)
1240 670 : CALL uppercase(int_exec_space_key)
1241 : params%use_gauxc_model = (TRIM(model_key) /= "" .AND. TRIM(model_key) /= "NONE" .AND. &
1242 670 : TRIM(model_key) /= TRIM(xc_fun_key))
1243 670 : use_skala_model = (INDEX(TRIM(model_key), "SKALA") > 0)
1244 670 : lwd_kernel_key = ADJUSTL(params%lwd_kernel)
1245 670 : CALL uppercase(lwd_kernel_key)
1246 670 : IF (TRIM(lwd_kernel_key) == "AUTO") THEN
1247 670 : params%lwd_kernel = "DEFAULT"
1248 : #if defined(GAUXC_HAS_CUTLASS)
1249 : IF (use_skala_model .AND. TRIM(params%int_exec_space) == "DEVICE") THEN
1250 : params%lwd_kernel = "SCHEME1-CUTLASS"
1251 : END IF
1252 : #endif
1253 : END IF
1254 670 : params%model_eval_name = model_name
1255 670 : IF (.NOT. params%use_gauxc_model) THEN
1256 : ! MODEL NONE and MODEL equal to FUNCTIONAL select conventional GauXC.
1257 510 : params%model_eval_name = "NONE"
1258 : END IF
1259 : IF (gapw_pseudopotentials .AND. params%use_gauxc_model .AND. .NOT. dft_control%qs_control%gapw_xc .AND. &
1260 670 : .NOT. gapw_paw_pseudopotentials .AND. para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
1261 2 : IF (scf_env%iter_count == 1) THEN
1262 : CALL cp_warn( &
1263 : __LOCATION__, &
1264 : "GauXC Skala with METHOD GAPW and GPW_TYPE pseudopotentials evaluates "// &
1265 : "the XC term directly on the molecular AO/valence density; no GAPW one-center "// &
1266 2 : "XC correction is used for those regular-grid kinds.")
1267 : END IF
1268 : END IF
1269 670 : IF (params%device_runtime_fill_fraction <= 0.0_dp .OR. params%device_runtime_fill_fraction > 1.0_dp) THEN
1270 : CALL cp_abort(__LOCATION__, &
1271 0 : "GAUXC%DEVICE_RUNTIME_FILL_FRACTION must be > 0 and <= 1.")
1272 : END IF
1273 670 : IF (atom_chunk_size < -1) THEN
1274 : CALL cp_abort(__LOCATION__, &
1275 0 : "GAUXC%MODEL_ATOM_CHUNK_SIZE must be -1, zero, or positive.")
1276 : END IF
1277 670 : IF (molecular_virial_debug) THEN
1278 48 : IF (molecular_virial_debug_dx <= 0.0_dp) THEN
1279 : CALL cp_abort(__LOCATION__, &
1280 0 : "GauXC MOLECULAR_VIRIAL_DEBUG_DX must be positive.")
1281 : END IF
1282 48 : molecular_virial = .TRUE.
1283 : END IF
1284 670 : need_xc_gradient = calculate_forces .OR. molecular_virial
1285 670 : IF (low_spin_scope) THEN
1286 96 : IF (.NOT. dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1287 0 : CPABORT("GauXC low spin ROKS requires METHOD GAPW.")
1288 : END IF
1289 96 : IF (.NOT. gauxc_gapw_all_all_electron(qs_kind_set)) THEN
1290 0 : CPABORT("GAPW low spin ROKS GauXC requires all-electron kinds.")
1291 : END IF
1292 96 : IF (is_periodic .OR. do_kpoints .OR. nimages /= 1) THEN
1293 0 : CPABORT("GAPW low spin ROKS GauXC requires an isolated Gamma-point calculation.")
1294 : END IF
1295 96 : IF (need_xc_gradient) THEN
1296 : CALL cp_abort(__LOCATION__, &
1297 0 : "GAPW low spin ROKS GauXC nuclear gradients and molecular virial are not implemented.")
1298 : END IF
1299 96 : IF (TRIM(int_exec_space_key) /= "HOST" .OR. TRIM(params%lb_exec_space) /= "HOST") THEN
1300 0 : CPABORT("GAPW low spin ROKS GauXC currently requires HOST execution.")
1301 : END IF
1302 96 : IF (params%use_gauxc_model .AND. .NOT. use_skala_model) THEN
1303 0 : CPABORT("GAPW low spin ROKS GauXC supports only the SKALA model or conventional XC.")
1304 : END IF
1305 96 : IF (use_skala_model .AND. TRIM(skala_runtime_key) == "MPI") THEN
1306 0 : CPABORT("GAPW low spin ROKS Skala requires SKALA_RUNTIME AUTO or SELF.")
1307 : END IF
1308 : END IF
1309 : CALL ensure_gauxc_periodic_reference_scope( &
1310 670 : dft_control, cell, qs_kind_set, do_kpoints, periodic_reference)
1311 670 : IF (is_periodic .AND. periodic_reference .AND. para_env%mepos == 0) THEN
1312 267 : IF (ASSOCIATED(scf_env)) THEN
1313 267 : IF (scf_env%iter_count == 1) THEN
1314 : CALL cp_warn( &
1315 : __LOCATION__, &
1316 : "GAUXC%PERIODIC_REFERENCE uses GauXC molecular quadrature for isolated validation "// &
1317 36 : "cells. Compact periodic materials require a dedicated periodic GauXC interface.")
1318 : END IF
1319 : END IF
1320 : END IF
1321 670 : IF (params%use_gauxc_model) THEN
1322 160 : IF (has_nlcc(qs_kind_set)) THEN
1323 : CALL cp_abort(__LOCATION__, &
1324 : "GauXC Skala with NLCC pseudopotentials is not implemented. "// &
1325 0 : "The frozen core density would need a SKALA-consistent feature definition.")
1326 : END IF
1327 : END IF
1328 670 : IF (params%use_gauxc_model) THEN
1329 : CALL set_gauxc_model_atom_chunk_env( &
1330 160 : atom_chunk_size, atom_chunk_size_explicit)
1331 160 : IF (.NOT. grid_explicit) params%grid_type = "SUPERFINE"
1332 160 : IF (.NOT. pruning_explicit) params%pruning_scheme = "UNPRUNED"
1333 :
1334 160 : grid_key = ADJUSTL(params%grid_type)
1335 160 : pruning_key = ADJUSTL(params%pruning_scheme)
1336 160 : CALL uppercase(grid_key)
1337 160 : CALL uppercase(pruning_key)
1338 160 : IF (use_skala_model .AND. need_xc_gradient .AND. &
1339 : (TRIM(grid_key) /= "SUPERFINE" .OR. TRIM(pruning_key) /= "UNPRUNED")) THEN
1340 : CALL cp_warn( &
1341 : __LOCATION__, &
1342 : "GauXC Skala nuclear gradients are sensitive to the GauXC molecular grid. "// &
1343 0 : "Use GRID SUPERFINE and PRUNING_SCHEME UNPRUNED for quantitative force checks.")
1344 : END IF
1345 160 : IF (TRIM(model_key) == "SKALA") THEN
1346 64 : model_name = ""
1347 64 : env_status = 1
1348 64 : IF (TRIM(int_exec_space_key) == "DEVICE") THEN
1349 0 : CALL GET_ENVIRONMENT_VARIABLE("SKALA_CUDA_MODEL", model_name, STATUS=env_status)
1350 : END IF
1351 64 : IF (env_status /= 0 .OR. LEN_TRIM(model_name) == 0) THEN
1352 64 : CALL GET_ENVIRONMENT_VARIABLE("SKALA_MODEL", model_name, STATUS=env_status)
1353 : END IF
1354 64 : IF (env_status /= 0 .OR. LEN_TRIM(model_name) == 0) THEN
1355 0 : IF (TRIM(int_exec_space_key) == "DEVICE") THEN
1356 : CALL cp_abort( &
1357 : __LOCATION__, &
1358 : "MODEL SKALA with DEVICE execution requires SKALA_CUDA_MODEL or "// &
1359 0 : "GAUXC_SKALA_MODEL")
1360 : ELSE
1361 0 : CPABORT("MODEL SKALA requires the GAUXC_SKALA_MODEL environment variable")
1362 : END IF
1363 : END IF
1364 64 : params%model_eval_name = model_name
1365 : END IF
1366 : END IF
1367 1244 : SELECT CASE (TRIM(skala_runtime_key))
1368 : CASE ("AUTO")
1369 594 : params%use_self_runtime = use_skala_model .AND. para_env%num_pe > 1 .AND. params%nspins > 1
1370 : CASE ("MPI")
1371 0 : params%use_self_runtime = .FALSE.
1372 : CASE ("SELF")
1373 96 : params%use_self_runtime = use_skala_model .AND. para_env%num_pe > 1
1374 : CASE DEFAULT
1375 670 : CALL cp_abort(__LOCATION__, "Unknown GAUXC%SKALA_RUNTIME value.")
1376 : END SELECT
1377 670 : IF (.NOT. use_skala_model) params%use_self_runtime = .FALSE.
1378 670 : params%gradient_runtime = gradient_runtime_key
1379 : IF (use_skala_model .AND. para_env%num_pe > 1 .AND. .NOT. params%use_self_runtime .AND. &
1380 670 : para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
1381 10 : IF (scf_env%iter_count == 1) THEN
1382 : CALL cp_warn( &
1383 : __LOCATION__, &
1384 : "GAUXC%SKALA_RUNTIME uses the MPI communicator for energy/VXC. "// &
1385 : "SKALA Torch atom chunks can be distributed across MPI ranks; "// &
1386 10 : "set GAUXC_ONEDFT_DISTRIBUTED_TORCH=0 to force rank-0 Torch inference.")
1387 : END IF
1388 : END IF
1389 :
1390 : ! Resolve the gradient route once the cache knows the actual basis angular momentum.
1391 670 : params%gapw_method = gapw_method
1392 670 : params%need_xc_gradient = need_xc_gradient
1393 :
1394 670 : IF (.NOT. ASSOCIATED(qs_env%gauxc_cache)) ALLOCATE (qs_env%gauxc_cache)
1395 670 : cache => qs_env%gauxc_cache
1396 : CALL gauxc_cache_init( &
1397 : cache, &
1398 : params, &
1399 : para_env, &
1400 : particle_set, &
1401 : qs_kind_set, &
1402 670 : gauxc_status)
1403 :
1404 670 : hdf5_output = (TRIM(output_path) /= "")
1405 670 : write_hdf5_output = hdf5_output .AND. para_env%mepos == 0 .AND. .NOT. external_density
1406 0 : IF (write_hdf5_output .AND. ASSOCIATED(scf_env)) THEN
1407 0 : write_hdf5_output = scf_env%iter_count == 1
1408 : END IF
1409 0 : IF (write_hdf5_output) THEN
1410 : CALL gauxc_write_molecule_hdf5( &
1411 : cache%molecule, &
1412 : output_path, &
1413 : "molecule.h5", &
1414 : "molecule", &
1415 0 : gauxc_status)
1416 0 : CALL gauxc_check_status(gauxc_status)
1417 : CALL gauxc_write_basisset_hdf5( &
1418 : cache%basisset, &
1419 : output_path, &
1420 : "basisset.h5", &
1421 : "basisset", &
1422 0 : gauxc_status)
1423 0 : CALL gauxc_check_status(gauxc_status)
1424 : END IF
1425 :
1426 670 : IF (qs_env%run_rtp) THEN
1427 0 : CPABORT("GAUXC XC energy currently does not support real-time propagation")
1428 : END IF
1429 :
1430 670 : exc_value = 0.0_dp
1431 :
1432 670 : IF (ASSOCIATED(matrix_vxc)) CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1433 670 : CALL dbcsr_allocate_matrix_set(matrix_vxc, params%nspins)
1434 :
1435 1340 : DO img = 1, nimages
1436 670 : IF (img > 1) THEN
1437 0 : CPABORT("UNIMPLEMENTED: Handling nimg>1 in k-point integration")
1438 : END IF
1439 670 : CALL dbcsr_to_dense(rho_ao(1, img), density_scalar, para_env)
1440 670 : CALL para_env%sum(density_scalar)
1441 670 : IF (params%nspins == 1) THEN
1442 : gauxc_xc_result = gauxc_compute_xc( &
1443 : cache%integrator, &
1444 : density_scalar, &
1445 : nspins=params%nspins, &
1446 : status=gauxc_status, &
1447 552 : model=TRIM(params%model_eval_name))
1448 552 : CALL gauxc_check_status(gauxc_status)
1449 552 : IF (need_xc_gradient) THEN
1450 56 : IF (params%use_fd_gradient) THEN
1451 : CALL gauxc_xc_gradient_fd( &
1452 : particle_set, qs_kind_set, density_scalar, params%nspins, params%model_eval_name, &
1453 : params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1454 : params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1455 : params%device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
1456 2 : exc_grad%exc_grad)
1457 54 : ELSE IF (params%use_gradient_self_runtime) THEN
1458 : exc_grad = gauxc_compute_xc_gradient( &
1459 : cache%gradient_integrator, &
1460 : density_scalar, &
1461 : nspins=params%nspins, &
1462 : natom=params%natom, &
1463 : status=gauxc_status, &
1464 48 : model=TRIM(params%model_eval_name))
1465 : ELSE
1466 : exc_grad = gauxc_compute_xc_gradient( &
1467 : cache%integrator, &
1468 : density_scalar, &
1469 : nspins=params%nspins, &
1470 : natom=params%natom, &
1471 : status=gauxc_status, &
1472 6 : model=TRIM(params%model_eval_name))
1473 : END IF
1474 56 : CALL gauxc_check_status(gauxc_status)
1475 56 : IF (calculate_forces) THEN
1476 : CALL add_gauxc_gradient_to_force( &
1477 : exc_grad%exc_grad, &
1478 : force, &
1479 : atomic_kind_set, &
1480 10 : para_env)
1481 : END IF
1482 56 : IF (molecular_virial) THEN
1483 48 : CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1484 : END IF
1485 56 : IF (molecular_virial_debug) THEN
1486 : CALL debug_gauxc_molecular_virial( &
1487 : exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, params%nspins, &
1488 : params%model_eval_name, params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1489 : params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1490 48 : params%device_runtime_fill_fraction, molecular_virial_debug_dx, para_env)
1491 : END IF
1492 56 : DEALLOCATE (exc_grad%exc_grad)
1493 : END IF
1494 : ELSE
1495 118 : CPASSERT(params%nspins == 2)
1496 : ! In here:
1497 : ! scalar <- rho_ao(1, :) + rho_ao(2, :)
1498 : ! zeta <- rho_ao(1, :) - rho_ao(2, :)
1499 118 : CALL dbcsr_to_dense(rho_ao(2, img), density_zeta, para_env)
1500 118 : CALL para_env%sum(density_zeta)
1501 : ! Do NOT reorder the following lines!
1502 83378 : density_scalar(:, :) = density_scalar(:, :) + density_zeta(:, :)
1503 : ! Factor two because the next line is evaluated after the above line.
1504 : ! We need to subtract density_zeta once to undo the above line and
1505 : ! a second time because that is what UKS requires.
1506 : ! This style lowers memory footprint.
1507 83378 : density_zeta(:, :) = density_scalar(:, :) - 2.0_dp*density_zeta(:, :)
1508 : gauxc_xc_result = gauxc_compute_xc( &
1509 : cache%integrator, &
1510 : density_scalar, &
1511 : density_zeta, &
1512 : params%nspins, &
1513 : gauxc_status, &
1514 118 : model=TRIM(params%model_eval_name))
1515 118 : CALL gauxc_check_status(gauxc_status)
1516 118 : IF (need_xc_gradient) THEN
1517 0 : IF (params%use_fd_gradient) THEN
1518 : CALL gauxc_xc_gradient_fd( &
1519 : particle_set, qs_kind_set, density_scalar, params%nspins, params%model_eval_name, &
1520 : params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1521 : params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1522 : params%device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
1523 0 : exc_grad%exc_grad, density_zeta=density_zeta)
1524 0 : ELSE IF (params%use_gradient_self_runtime) THEN
1525 : exc_grad = gauxc_compute_xc_gradient( &
1526 : cache%gradient_integrator, &
1527 : density_scalar, &
1528 : density_zeta, &
1529 : params%nspins, &
1530 : params%natom, &
1531 : gauxc_status, &
1532 0 : model=TRIM(params%model_eval_name))
1533 : ELSE
1534 : exc_grad = gauxc_compute_xc_gradient( &
1535 : cache%integrator, &
1536 : density_scalar, &
1537 : density_zeta, &
1538 : params%nspins, &
1539 : params%natom, &
1540 : gauxc_status, &
1541 0 : model=TRIM(params%model_eval_name))
1542 : END IF
1543 0 : CALL gauxc_check_status(gauxc_status)
1544 0 : IF (calculate_forces) THEN
1545 : CALL add_gauxc_gradient_to_force( &
1546 : exc_grad%exc_grad, &
1547 : force, &
1548 : atomic_kind_set, &
1549 0 : para_env)
1550 : END IF
1551 0 : IF (molecular_virial) THEN
1552 0 : CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1553 : END IF
1554 0 : IF (molecular_virial_debug) THEN
1555 : CALL debug_gauxc_molecular_virial( &
1556 : exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, params%nspins, &
1557 : params%model_eval_name, params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1558 : params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1559 : params%device_runtime_fill_fraction, molecular_virial_debug_dx, para_env, &
1560 0 : density_zeta=density_zeta)
1561 : END IF
1562 0 : DEALLOCATE (exc_grad%exc_grad)
1563 : END IF
1564 : END IF
1565 :
1566 670 : exc_value = exc_value + gauxc_xc_result%exc
1567 :
1568 1340 : IF (params%nspins == 1) THEN
1569 552 : IF (img == 1) THEN
1570 552 : matrix_vxc(1) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(1, img))
1571 : ELSE
1572 0 : CPABORT("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1573 : END IF
1574 : ELSE
1575 118 : CPASSERT(params%nspins == 2)
1576 : ! Transform derivatives from total/spin density back to alpha/beta channels.
1577 118 : vxc_zeta_tmp = dense_to_dbcsr(gauxc_xc_result%vxc_zeta, rho_ao(1, img))
1578 118 : IF (img == 1) THEN
1579 354 : DO ispin = 1, 2
1580 236 : matrix_vxc(ispin) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(ispin, 1))
1581 : CALL dbcsr_add( &
1582 : matrix_vxc(ispin)%matrix, &
1583 : vxc_zeta_tmp%matrix, &
1584 : 1.0_dp, &
1585 : ! 1.0 for ispin==1, -1.0 for ispin==2
1586 354 : 1.0_dp - REAL(ispin - 1, dp)*2.0_dp)
1587 : END DO
1588 : ELSE
1589 0 : CPABORT("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1590 : END IF
1591 118 : CALL dbcsr_release(vxc_zeta_tmp%matrix)
1592 118 : DEALLOCATE (vxc_zeta_tmp%matrix)
1593 : END IF
1594 : END DO
1595 :
1596 670 : DEALLOCATE (density_scalar)
1597 670 : IF (ALLOCATED(density_zeta)) DEALLOCATE (density_zeta)
1598 670 : DEALLOCATE (gauxc_xc_result%vxc_scalar)
1599 670 : IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
1600 :
1601 670 : IF (external_density) THEN
1602 64 : exc_external = exc_value
1603 64 : matrix_vxc_external => matrix_vxc
1604 : ELSE
1605 606 : energy%exc = exc_value
1606 606 : CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
1607 : END IF
1608 1458 : DO ispin = 1, params%nspins
1609 1458 : CALL dbcsr_finalize(matrix_vxc(ispin)%matrix)
1610 : END DO
1611 :
1612 1340 : END SUBROUTINE apply_gauxc
1613 :
1614 : END MODULE xc_gauxc_functional
|
