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