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 :: gpw_type_forced, 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 : gpw_type_forced=gpw_type_forced, &
401 : paw_atom=paw_atom, &
402 12 : sgp_potential=sgp_potential)
403 : IF ((ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) .AND. &
404 22 : paw_atom .AND. .NOT. gpw_type_forced) THEN
405 : has_paw_pseudopotentials = .TRUE.
406 : EXIT
407 : END IF
408 : END DO
409 :
410 10 : END FUNCTION gauxc_gapw_has_paw_pseudopotentials
411 :
412 : ! **************************************************************************************************
413 : !> \brief Check the current periodic scope of the CP2K-GauXC bridge
414 : !> \param dft_control ...
415 : !> \param cell ...
416 : !> \param qs_kind_set ...
417 : !> \param do_kpoints ...
418 : !> \param periodic_reference ...
419 : !> \note This path keeps isolated validation cells usable under PERIODIC XYZ.
420 : !> It intentionally does not implement compact periodic GauXC quadrature.
421 : ! **************************************************************************************************
422 378 : SUBROUTINE ensure_gauxc_periodic_reference_scope( &
423 : dft_control, cell, qs_kind_set, do_kpoints, periodic_reference)
424 : TYPE(dft_control_type), POINTER :: dft_control
425 : TYPE(cell_type), POINTER :: cell
426 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
427 : LOGICAL, INTENT(IN) :: do_kpoints, periodic_reference
428 :
429 : INTEGER :: ikind
430 : LOGICAL :: is_periodic
431 : TYPE(gth_potential_type), POINTER :: gth_potential
432 : TYPE(sgp_potential_type), POINTER :: sgp_potential
433 :
434 378 : CPASSERT(ASSOCIATED(dft_control))
435 378 : CPASSERT(ASSOCIATED(qs_kind_set))
436 :
437 378 : is_periodic = .FALSE.
438 414 : IF (ASSOCIATED(cell)) is_periodic = ANY(cell%perd /= 0)
439 :
440 378 : IF (do_kpoints) THEN
441 : CALL cp_abort(__LOCATION__, &
442 : "GauXC currently supports only Gamma-only density matrices in CP2K. "// &
443 0 : "Periodic k-point density matrices require a dedicated GauXC periodic interface.")
444 : END IF
445 378 : IF (dft_control%nimages /= 1) THEN
446 : CALL cp_abort(__LOCATION__, &
447 : "GauXC currently supports only a single AO image in CP2K. "// &
448 0 : "Periodic neighbour-cell AO blocks require a dedicated GauXC periodic interface.")
449 : END IF
450 378 : IF (.NOT. is_periodic) RETURN
451 :
452 366 : IF (.NOT. periodic_reference) THEN
453 : CALL cp_abort(__LOCATION__, &
454 : "Periodic GauXC calculations in CP2K require GAUXC%PERIODIC_REFERENCE T. "// &
455 : "This opt-in documents that the current path is only an isolated-cell, "// &
456 : "Gamma-only, single-image METHOD GPW reference path using GauXC molecular "// &
457 0 : "quadrature, not a dedicated periodic GauXC interface.")
458 : END IF
459 :
460 1464 : IF (.NOT. ALL(cell%perd == 1)) THEN
461 : CALL cp_abort(__LOCATION__, &
462 : "The current GauXC isolated-cell reference path supports only PERIODIC XYZ. "// &
463 0 : "Partial periodicity requires a dedicated GauXC periodic interface.")
464 : END IF
465 366 : IF (.NOT. dft_control%qs_control%gpw) THEN
466 : CALL cp_abort(__LOCATION__, &
467 : "The current GauXC isolated-cell reference path is limited to METHOD GPW with GTH "// &
468 0 : "pseudopotentials. GAPW, GAPW_XC, and other QS methods are not supported here.")
469 : END IF
470 :
471 846 : DO ikind = 1, SIZE(qs_kind_set)
472 480 : NULLIFY (gth_potential, sgp_potential)
473 : CALL get_qs_kind(qs_kind_set(ikind), &
474 : gth_potential=gth_potential, &
475 480 : sgp_potential=sgp_potential)
476 846 : IF (.NOT. ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
477 : CALL cp_abort(__LOCATION__, &
478 : "The current GauXC isolated-cell reference path is limited to GTH pseudopotentials. "// &
479 0 : "Use non-periodic all-electron GAPW validation for molecular GAPW cases.")
480 : END IF
481 : END DO
482 :
483 : END SUBROUTINE ensure_gauxc_periodic_reference_scope
484 :
485 : ! **************************************************************************************************
486 : !> \brief adds a replicated GauXC energy gradient to the local CP2K force accumulator
487 : !> \param exc_grad ...
488 : !> \param force ...
489 : !> \param atomic_kind_set ...
490 : !> \param para_env ...
491 : ! **************************************************************************************************
492 4 : SUBROUTINE add_gauxc_gradient_to_force(exc_grad, force, atomic_kind_set, para_env)
493 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
494 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
495 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
496 : TYPE(mp_para_env_type), POINTER :: para_env
497 :
498 : INTEGER :: ia, iatom, ikind, natom_kind
499 : TYPE(atomic_kind_type), POINTER :: atomic_kind
500 :
501 4 : CPASSERT(ASSOCIATED(force))
502 4 : CPASSERT(ASSOCIATED(atomic_kind_set))
503 :
504 4 : IF (para_env%mepos /= 0) RETURN
505 :
506 5 : DO ikind = 1, SIZE(atomic_kind_set, 1)
507 3 : atomic_kind => atomic_kind_set(ikind)
508 3 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
509 11 : DO ia = 1, natom_kind
510 6 : iatom = atomic_kind%atom_list(ia)
511 : force(ikind)%rho_elec(:, ia) = force(ikind)%rho_elec(:, ia) + &
512 27 : exc_grad(3*iatom - 2:3*iatom)
513 : END DO
514 : END DO
515 :
516 : END SUBROUTINE add_gauxc_gradient_to_force
517 :
518 : ! **************************************************************************************************
519 : !> \brief compute a GauXC XC energy for diagnostic finite differences
520 : !> \param particle_set_eval ...
521 : !> \param qs_kind_set ...
522 : !> \param density_scalar ...
523 : !> \param nspins ...
524 : !> \param model_name ...
525 : !> \param xc_fun_name ...
526 : !> \param grid_type ...
527 : !> \param radial_quadrature ...
528 : !> \param pruning_scheme ...
529 : !> \param lb_exec_space ...
530 : !> \param int_exec_space ...
531 : !> \param lwd_kernel ...
532 : !> \param batch_size ...
533 : !> \param device_runtime_fill_fraction ...
534 : !> \param exc ...
535 : !> \param density_zeta ...
536 : ! **************************************************************************************************
537 0 : SUBROUTINE gauxc_xc_energy_for_particles( &
538 0 : particle_set_eval, qs_kind_set, density_scalar, nspins, model_name, &
539 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
540 0 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, exc, density_zeta)
541 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_eval
542 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
543 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
544 : INTEGER, INTENT(IN) :: nspins
545 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
546 : pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
547 : INTEGER, INTENT(IN) :: batch_size
548 : REAL(KIND=dp), INTENT(IN) :: device_runtime_fill_fraction
549 : REAL(KIND=dp), INTENT(OUT) :: exc
550 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
551 : OPTIONAL :: density_zeta
552 :
553 : TYPE(cp_gauxc_basisset_type) :: gauxc_basis_fd
554 : TYPE(cp_gauxc_grid_type) :: gauxc_grid_fd
555 : TYPE(cp_gauxc_integrator_type) :: gauxc_integrator_fd
556 : TYPE(cp_gauxc_molecule_type) :: gauxc_mol_fd
557 : TYPE(cp_gauxc_status_type) :: gauxc_status
558 0 : TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
559 :
560 0 : gauxc_mol_fd = gauxc_create_molecule(particle_set_eval, gauxc_status)
561 0 : CALL gauxc_check_status(gauxc_status)
562 0 : gauxc_basis_fd = gauxc_create_basisset(qs_kind_set, particle_set_eval, gauxc_status)
563 0 : CALL gauxc_check_status(gauxc_status)
564 : gauxc_grid_fd = gauxc_create_grid( &
565 : gauxc_mol_fd, &
566 : gauxc_basis_fd, &
567 : grid_type, &
568 : radial_quadrature, &
569 : pruning_scheme, &
570 : lb_exec_space, &
571 : batch_size, &
572 : device_runtime_fill_fraction, &
573 : gauxc_status, &
574 0 : mpi_comm=mp_comm_self%get_handle())
575 0 : CALL gauxc_check_status(gauxc_status)
576 : gauxc_integrator_fd = gauxc_create_integrator( &
577 : TRIM(xc_fun_name), &
578 : gauxc_grid_fd, &
579 : int_exec_space, &
580 : lwd_kernel, &
581 : nspins, &
582 0 : gauxc_status)
583 0 : CALL gauxc_check_status(gauxc_status)
584 :
585 0 : IF (nspins == 1) THEN
586 : gauxc_xc_result = gauxc_compute_xc( &
587 : gauxc_integrator_fd, &
588 : density_scalar, &
589 : nspins=nspins, &
590 : status=gauxc_status, &
591 0 : model=TRIM(model_name))
592 : ELSE
593 0 : CPASSERT(nspins == 2)
594 0 : CPASSERT(PRESENT(density_zeta))
595 : gauxc_xc_result = gauxc_compute_xc( &
596 : gauxc_integrator_fd, &
597 : density_scalar, &
598 : density_zeta, &
599 : nspins, &
600 : gauxc_status, &
601 0 : model=TRIM(model_name))
602 : END IF
603 0 : CALL gauxc_check_status(gauxc_status)
604 0 : exc = gauxc_xc_result%exc
605 :
606 0 : IF (ALLOCATED(gauxc_xc_result%vxc_scalar)) DEALLOCATE (gauxc_xc_result%vxc_scalar)
607 0 : IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
608 :
609 0 : CALL gauxc_destroy_integrator(gauxc_integrator_fd, gauxc_status)
610 0 : CALL gauxc_check_status(gauxc_status)
611 0 : CALL gauxc_destroy_grid(gauxc_grid_fd, gauxc_status)
612 0 : CALL gauxc_check_status(gauxc_status)
613 0 : CALL gauxc_destroy_basisset(gauxc_basis_fd, gauxc_status)
614 0 : CALL gauxc_check_status(gauxc_status)
615 0 : CALL gauxc_destroy_molecule(gauxc_mol_fd, gauxc_status)
616 0 : CALL gauxc_check_status(gauxc_status)
617 :
618 0 : END SUBROUTINE gauxc_xc_energy_for_particles
619 :
620 : ! **************************************************************************************************
621 : !> \brief compute a finite-difference GauXC XC nuclear gradient at fixed density
622 : !> \param particle_set ...
623 : !> \param qs_kind_set ...
624 : !> \param density_scalar ...
625 : !> \param nspins ...
626 : !> \param model_name ...
627 : !> \param xc_fun_name ...
628 : !> \param grid_type ...
629 : !> \param radial_quadrature ...
630 : !> \param pruning_scheme ...
631 : !> \param lb_exec_space ...
632 : !> \param int_exec_space ...
633 : !> \param lwd_kernel ...
634 : !> \param batch_size ...
635 : !> \param device_runtime_fill_fraction ...
636 : !> \param dx ...
637 : !> \param para_env ...
638 : !> \param exc_grad ...
639 : !> \param density_zeta ...
640 : ! **************************************************************************************************
641 0 : SUBROUTINE gauxc_xc_gradient_fd( &
642 0 : particle_set, qs_kind_set, density_scalar, nspins, model_name, &
643 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
644 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, dx, para_env, exc_grad, &
645 0 : density_zeta)
646 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
647 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
648 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
649 : INTEGER, INTENT(IN) :: nspins
650 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
651 : pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
652 : INTEGER, INTENT(IN) :: batch_size
653 : REAL(KIND=dp), INTENT(IN) :: device_runtime_fill_fraction, dx
654 : TYPE(mp_para_env_type), POINTER :: para_env
655 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
656 : INTENT(OUT) :: exc_grad
657 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
658 : OPTIONAL :: density_zeta
659 :
660 : INTEGER :: iatom, idir
661 : REAL(KIND=dp) :: xc_minus, xc_plus
662 0 : TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus
663 :
664 0 : CPASSERT(ASSOCIATED(particle_set))
665 0 : CPASSERT(dx > 0.0_dp)
666 :
667 0 : ALLOCATE (exc_grad(3*SIZE(particle_set)))
668 0 : exc_grad = 0.0_dp
669 :
670 0 : IF (para_env%mepos == 0) THEN
671 0 : ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set)))
672 :
673 0 : DO iatom = 1, SIZE(particle_set)
674 0 : DO idir = 1, 3
675 0 : particle_set_minus = particle_set
676 0 : particle_set_plus = particle_set
677 0 : particle_set_minus(iatom)%r(idir) = particle_set_minus(iatom)%r(idir) - dx
678 0 : particle_set_plus(iatom)%r(idir) = particle_set_plus(iatom)%r(idir) + dx
679 0 : IF (PRESENT(density_zeta)) THEN
680 : CALL gauxc_xc_energy_for_particles( &
681 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
682 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
683 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus, &
684 0 : density_zeta=density_zeta)
685 : CALL gauxc_xc_energy_for_particles( &
686 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
687 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
688 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus, &
689 0 : density_zeta=density_zeta)
690 : ELSE
691 : CALL gauxc_xc_energy_for_particles( &
692 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
693 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
694 0 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus)
695 : CALL gauxc_xc_energy_for_particles( &
696 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
697 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
698 0 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus)
699 : END IF
700 0 : exc_grad(3*iatom - 3 + idir) = (xc_plus - xc_minus)/(2.0_dp*dx)
701 : END DO
702 : END DO
703 :
704 0 : DEALLOCATE (particle_set_minus, particle_set_plus)
705 : END IF
706 :
707 0 : CALL para_env%bcast(exc_grad, 0)
708 :
709 0 : END SUBROUTINE gauxc_xc_gradient_fd
710 :
711 : ! **************************************************************************************************
712 : !> \brief finite-difference check of the molecular GauXC XC virial diagnostic
713 : !> \param exc_grad ...
714 : !> \param particle_set ...
715 : !> \param qs_kind_set ...
716 : !> \param density_scalar ...
717 : !> \param nspins ...
718 : !> \param model_name ...
719 : !> \param xc_fun_name ...
720 : !> \param grid_type ...
721 : !> \param radial_quadrature ...
722 : !> \param pruning_scheme ...
723 : !> \param lb_exec_space ...
724 : !> \param int_exec_space ...
725 : !> \param lwd_kernel ...
726 : !> \param batch_size ...
727 : !> \param device_runtime_fill_fraction ...
728 : !> \param dx ...
729 : !> \param para_env ...
730 : !> \param density_zeta ...
731 : ! **************************************************************************************************
732 0 : SUBROUTINE debug_gauxc_molecular_virial( &
733 0 : exc_grad, particle_set, qs_kind_set, density_scalar, nspins, model_name, &
734 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
735 0 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, dx, para_env, density_zeta)
736 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
737 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
738 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
739 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
740 : INTEGER, INTENT(IN) :: nspins
741 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
742 : pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
743 : INTEGER, INTENT(IN) :: batch_size
744 : REAL(KIND=dp), INTENT(IN) :: device_runtime_fill_fraction, dx
745 : TYPE(mp_para_env_type), POINTER :: para_env
746 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
747 : OPTIONAL :: density_zeta
748 :
749 : INTEGER :: iatom, iw
750 : REAL(KIND=dp) :: analytic_trace, diff_trace, &
751 : numerical_trace, xc_minus, xc_plus
752 : REAL(KIND=dp), DIMENSION(3) :: center, displacement, grad
753 0 : TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus
754 :
755 0 : CPASSERT(ASSOCIATED(particle_set))
756 0 : CPASSERT(SIZE(exc_grad) == 3*SIZE(particle_set))
757 :
758 0 : IF (para_env%mepos /= 0) RETURN
759 :
760 0 : center = 0.0_dp
761 0 : DO iatom = 1, SIZE(particle_set)
762 0 : center = center + particle_set(iatom)%r
763 : END DO
764 0 : center = center/REAL(SIZE(particle_set), dp)
765 :
766 0 : ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set)))
767 0 : particle_set_minus = particle_set
768 0 : particle_set_plus = particle_set
769 :
770 0 : analytic_trace = 0.0_dp
771 0 : DO iatom = 1, SIZE(particle_set)
772 0 : grad = exc_grad(3*iatom - 2:3*iatom)
773 0 : displacement = particle_set(iatom)%r - center
774 0 : analytic_trace = analytic_trace + DOT_PRODUCT(grad, displacement)
775 0 : particle_set_minus(iatom)%r = center + (1.0_dp - dx)*displacement
776 0 : particle_set_plus(iatom)%r = center + (1.0_dp + dx)*displacement
777 : END DO
778 0 : analytic_trace = analytic_trace/3.0_dp
779 :
780 0 : IF (PRESENT(density_zeta)) THEN
781 : CALL gauxc_xc_energy_for_particles( &
782 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
783 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
784 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus, &
785 0 : density_zeta=density_zeta)
786 : CALL gauxc_xc_energy_for_particles( &
787 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
788 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
789 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus, &
790 0 : density_zeta=density_zeta)
791 : ELSE
792 : CALL gauxc_xc_energy_for_particles( &
793 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
794 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
795 0 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus)
796 : CALL gauxc_xc_energy_for_particles( &
797 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
798 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
799 0 : int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus)
800 : END IF
801 :
802 0 : numerical_trace = (xc_plus - xc_minus)/(2.0_dp*dx)/3.0_dp
803 0 : diff_trace = analytic_trace - numerical_trace
804 :
805 0 : iw = cp_logger_get_default_io_unit()
806 0 : IF (iw > 0) THEN
807 : WRITE (UNIT=iw, FMT="(/,T2,A,1X,ES11.4)") &
808 0 : "GAUXC| Molecular XC virial finite-difference dx", dx
809 : WRITE (UNIT=iw, FMT="(T2,A,3(1X,ES19.11))") &
810 0 : "GAUXC| Molecular XC virial FD 1/3 Trace", &
811 0 : analytic_trace, numerical_trace, diff_trace
812 : END IF
813 :
814 0 : DEALLOCATE (particle_set_minus, particle_set_plus)
815 :
816 : END SUBROUTINE debug_gauxc_molecular_virial
817 :
818 : ! **************************************************************************************************
819 : !> \brief prints a force-based molecular XC virial diagnostic from GauXC gradients
820 : !> \param exc_grad ...
821 : !> \param particle_set ...
822 : !> \param para_env ...
823 : ! **************************************************************************************************
824 0 : SUBROUTINE print_gauxc_molecular_virial(exc_grad, particle_set, para_env)
825 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
826 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
827 : TYPE(mp_para_env_type), POINTER :: para_env
828 :
829 : CHARACTER(len=1), DIMENSION(3), PARAMETER :: label = ["x", "y", "z"]
830 :
831 : INTEGER :: i, iatom, iw, j
832 : REAL(KIND=dp), DIMENSION(3) :: center, displacement, grad, grad_sum
833 : REAL(KIND=dp), DIMENSION(3, 3) :: molecular_virial
834 :
835 0 : CPASSERT(ASSOCIATED(particle_set))
836 0 : CPASSERT(SIZE(exc_grad) == 3*SIZE(particle_set))
837 :
838 0 : IF (para_env%mepos /= 0) RETURN
839 :
840 0 : center = 0.0_dp
841 0 : DO iatom = 1, SIZE(particle_set)
842 0 : center = center + particle_set(iatom)%r
843 : END DO
844 0 : center = center/REAL(SIZE(particle_set), dp)
845 :
846 0 : grad_sum = 0.0_dp
847 0 : molecular_virial = 0.0_dp
848 0 : DO iatom = 1, SIZE(particle_set)
849 0 : grad = exc_grad(3*iatom - 2:3*iatom)
850 0 : displacement = particle_set(iatom)%r - center
851 0 : grad_sum = grad_sum + grad
852 0 : DO i = 1, 3
853 0 : DO j = 1, 3
854 0 : molecular_virial(i, j) = molecular_virial(i, j) + grad(i)*displacement(j)
855 : END DO
856 : END DO
857 : END DO
858 :
859 0 : iw = cp_logger_get_default_io_unit()
860 0 : IF (iw <= 0) RETURN
861 :
862 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
863 0 : "GAUXC| Molecular XC gradient virial diagnostic [a.u.]"
864 0 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,T40,A,T60,A)") "GAUXC|", "x", "y", "z"
865 0 : DO i = 1, 3
866 : WRITE (UNIT=iw, FMT="(T2,A,1X,A1,3(1X,ES19.11))") &
867 0 : "GAUXC|", label(i), molecular_virial(i, :)
868 : END DO
869 : WRITE (UNIT=iw, FMT="(T2,A,1X,ES19.11)") &
870 0 : "GAUXC| Molecular XC gradient virial 1/3 Trace", &
871 0 : (molecular_virial(1, 1) + molecular_virial(2, 2) + molecular_virial(3, 3))/3.0_dp
872 : WRITE (UNIT=iw, FMT="(T2,A,3(1X,ES19.11))") &
873 0 : "GAUXC| Molecular XC gradient sum", grad_sum
874 : WRITE (UNIT=iw, FMT="(T2,A)") &
875 0 : "GAUXC| Diagnostic only; this is not an analytical periodic stress tensor."
876 :
877 : END SUBROUTINE print_gauxc_molecular_virial
878 :
879 : ! **************************************************************************************************
880 : !> \brief Return information about the Skala functional
881 : !> \param functional section containing the SKALA subsection
882 : !> \param lsd if you are using lsd or lda
883 : !> \param reference the reference to the article where the functional is explained
884 : !> \param shortform the short definition of the functional
885 : !> \param needs the flags corresponding to the inputs needed by this
886 : !> functional are set to true (the flags not needed aren't touched)
887 : !> \param max_deriv the maximal derivative available
888 : ! **************************************************************************************************
889 352 : SUBROUTINE skala_info(functional, lsd, reference, shortform, needs, max_deriv)
890 : TYPE(section_vals_type), POINTER :: functional
891 : LOGICAL, INTENT(in) :: lsd
892 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
893 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
894 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
895 :
896 : CHARACTER(len=default_path_length) :: model_key, model_name
897 : CHARACTER(len=default_string_length) :: xc_fun_name
898 : LOGICAL :: native_grid
899 :
900 176 : CALL section_vals_val_get(functional, "FUNCTIONAL", c_val=xc_fun_name)
901 176 : CALL section_vals_val_get(functional, "MODEL", c_val=model_name)
902 176 : CALL section_vals_val_get(functional, "NATIVE_GRID", l_val=native_grid)
903 176 : model_key = ADJUSTL(model_name)
904 176 : CALL uppercase(model_key)
905 :
906 176 : IF (PRESENT(reference)) THEN
907 8 : IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "") THEN
908 1 : reference = "Functional computed by GauXC (underlying: "//TRIM(xc_fun_name)//")"
909 : ELSE
910 7 : reference = "Functional computed by GauXC OneDFT model "//TRIM(model_name)
911 : END IF
912 : END IF
913 176 : IF (PRESENT(shortform)) THEN
914 8 : IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "") THEN
915 1 : shortform = "GAUXC ("//TRIM(xc_fun_name)//")"
916 : ELSE
917 7 : shortform = "GAUXC OneDFT"
918 : END IF
919 : END IF
920 176 : IF (PRESENT(needs)) THEN
921 168 : IF (native_grid .AND. TRIM(model_key) /= "NONE" .AND. TRIM(model_key) /= "") THEN
922 96 : IF (lsd) THEN
923 16 : needs%rho_spin = .TRUE.
924 16 : needs%drho_spin = .TRUE.
925 16 : needs%tau_spin = .TRUE.
926 : ELSE
927 80 : needs%rho = .TRUE.
928 80 : needs%drho = .TRUE.
929 80 : needs%tau = .TRUE.
930 : END IF
931 : ELSE
932 72 : needs%rho = .TRUE.
933 72 : IF (lsd) THEN
934 12 : needs%rho_spin = .TRUE.
935 : END IF
936 : END IF
937 : END IF
938 176 : IF (PRESENT(max_deriv)) max_deriv = 1
939 :
940 176 : END SUBROUTINE skala_info
941 :
942 : ! GauXC uses replicated dense density and VXC matrices. The DBCSR density matrix
943 : ! is distributed over MPI ranks, so apply_gauxc allreduces the dense copy before
944 : ! passing it to GauXC.
945 :
946 : ! **************************************************************************************************
947 : !> \brief ...
948 : !> \param qs_env ...
949 : !> \param xc_section ...
950 : !> \param calculate_forces ...
951 : ! **************************************************************************************************
952 378 : SUBROUTINE apply_gauxc(qs_env, xc_section, calculate_forces)
953 : TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
954 : TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
955 : LOGICAL, INTENT(IN) :: calculate_forces
956 :
957 : CHARACTER(len=*), PARAMETER :: gapw_xc_abort_message = &
958 : "GauXC with METHOD GAPW_XC is not supported yet. "// &
959 : "The GAPW_XC one-center XC correction needs a dedicated GauXC design.", &
960 : nonlocal_vdw_abort_message = &
961 : "GauXC does not support non-local VDW_POTENTIAL corrections. "// &
962 : "Use an additive PAIR_POTENTIAL dispersion correction or disable GauXC."
963 : REAL(KIND=dp), PARAMETER :: gapw_fd_gradient_dx = 1.0E-4_dp
964 :
965 : CHARACTER(len=default_path_length) :: model_key, model_name, output_path
966 : CHARACTER(len=default_string_length) :: grid_key, grid_type, int_exec_space, lb_exec_space, &
967 : lwd_kernel, onedft_gradient_runtime, onedft_gradient_runtime_key, pruning_key, &
968 : pruning_scheme, radial_quadrature, skala_runtime, skala_runtime_key, xc_fun_name
969 : INTEGER :: batch_size, env_status, img, ispin, &
970 : natom, nimages, nspins, &
971 : onedft_atom_chunk_size
972 : LOGICAL :: do_kpoints, gapw_paw_pseudopotentials, gapw_pseudopotentials, grid_explicit, &
973 : hdf5_output, is_periodic, molecular_virial, molecular_virial_debug, &
974 : onedft_atom_chunk_size_explicit, periodic_reference, pruning_explicit, use_fd_gradient, &
975 : use_gradient_mpi_runtime, use_gradient_self_runtime, use_onedft, use_self_runtime, &
976 : use_skala_model, write_hdf5_output
977 : REAL(KIND=dp) :: device_runtime_fill_fraction, &
978 : molecular_virial_debug_dx
979 378 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: density_scalar, density_zeta
980 378 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
981 : TYPE(cell_type), POINTER :: cell
982 : TYPE(cp_gauxc_basisset_type) :: gauxc_basis
983 : TYPE(cp_gauxc_grid_type) :: gauxc_gradient_grid_result, &
984 : gauxc_grid_result
985 : TYPE(cp_gauxc_integrator_type) :: gauxc_gradient_integrator_result, gauxc_integrator_result
986 : TYPE(cp_gauxc_molecule_type) :: gauxc_mol
987 : TYPE(cp_gauxc_status_type) :: gauxc_status
988 378 : TYPE(cp_gauxc_xc_gradient_type) :: exc_grad
989 378 : TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
990 : TYPE(dbcsr_p_type) :: vxc_zeta_tmp
991 378 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
992 378 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
993 : TYPE(dft_control_type), POINTER :: dft_control
994 : TYPE(mp_para_env_type), POINTER :: para_env
995 378 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
996 : TYPE(qs_energy_type), POINTER :: energy
997 378 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
998 378 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
999 : TYPE(qs_ks_env_type), POINTER :: ks_env
1000 : TYPE(qs_rho_type), POINTER :: rho, rho_use, rho_xc
1001 : TYPE(qs_scf_env_type), POINTER :: scf_env
1002 : TYPE(section_vals_type), POINTER :: gauxc_functional_section
1003 :
1004 : NULLIFY ( &
1005 : atomic_kind_set, &
1006 378 : cell, &
1007 378 : dft_control, &
1008 378 : energy, &
1009 378 : force, &
1010 378 : ks_env, &
1011 378 : matrix_vxc, &
1012 378 : para_env, &
1013 378 : particle_set, &
1014 378 : qs_kind_set, &
1015 378 : rho, &
1016 378 : rho_use, &
1017 378 : rho_xc, &
1018 378 : rho_ao, &
1019 378 : scf_env)
1020 :
1021 : CALL get_qs_env( &
1022 : qs_env, &
1023 : cell=cell, &
1024 : dft_control=dft_control, &
1025 : do_kpoints=do_kpoints, &
1026 : energy=energy, &
1027 : ks_env=ks_env, &
1028 : matrix_vxc=matrix_vxc, &
1029 : natom=natom, &
1030 : atomic_kind_set=atomic_kind_set, &
1031 : force=force, &
1032 : para_env=para_env, &
1033 : particle_set=particle_set, &
1034 : qs_kind_set=qs_kind_set, &
1035 : rho=rho, &
1036 : rho_xc=rho_xc, &
1037 378 : scf_env=scf_env)
1038 :
1039 378 : IF (dft_control%qs_control%gapw_xc) THEN
1040 0 : CPABORT(gapw_xc_abort_message)
1041 : END IF
1042 : gapw_pseudopotentials = dft_control%qs_control%gapw .AND. &
1043 378 : gauxc_gapw_has_pseudopotentials(qs_kind_set)
1044 : gapw_paw_pseudopotentials = dft_control%qs_control%gapw .AND. &
1045 378 : gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
1046 378 : CPASSERT(ASSOCIATED(rho))
1047 378 : rho_use => rho
1048 : CALL qs_rho_get( &
1049 : rho_use, &
1050 378 : rho_ao_kp=rho_ao)
1051 :
1052 378 : nimages = dft_control%nimages
1053 378 : nspins = dft_control%nspins
1054 378 : is_periodic = .FALSE.
1055 414 : IF (ASSOCIATED(cell)) is_periodic = ANY(cell%perd /= 0)
1056 :
1057 378 : IF (ASSOCIATED(qs_env%dispersion_env)) THEN
1058 378 : IF (qs_env%dispersion_env%type == xc_vdw_fun_nonloc) THEN
1059 0 : CPABORT(nonlocal_vdw_abort_message)
1060 : END IF
1061 : END IF
1062 378 : NULLIFY (vxc_zeta_tmp%matrix)
1063 :
1064 378 : gauxc_functional_section => get_gauxc_functional(xc_section)
1065 : CALL section_vals_val_get( &
1066 : gauxc_functional_section, &
1067 : "FUNCTIONAL", &
1068 378 : c_val=xc_fun_name)
1069 : CALL section_vals_val_get( &
1070 : gauxc_functional_section, &
1071 : "MODEL", &
1072 378 : c_val=model_name)
1073 : CALL section_vals_val_get( &
1074 : gauxc_functional_section, &
1075 : "GRID", &
1076 : c_val=grid_type, &
1077 378 : explicit=grid_explicit)
1078 : CALL section_vals_val_get( &
1079 : gauxc_functional_section, &
1080 : "RADIAL_QUADRATURE", &
1081 378 : c_val=radial_quadrature)
1082 : CALL section_vals_val_get( &
1083 : gauxc_functional_section, &
1084 : "PRUNING_SCHEME", &
1085 : c_val=pruning_scheme, &
1086 378 : explicit=pruning_explicit)
1087 : CALL section_vals_val_get( &
1088 : gauxc_functional_section, &
1089 : "BATCH_SIZE", &
1090 378 : i_val=batch_size)
1091 : CALL section_vals_val_get( &
1092 : gauxc_functional_section, &
1093 : "DEVICE_RUNTIME_FILL_FRACTION", &
1094 378 : r_val=device_runtime_fill_fraction)
1095 : CALL section_vals_val_get( &
1096 : gauxc_functional_section, &
1097 : "ONEDFT_ATOM_CHUNK_SIZE", &
1098 : i_val=onedft_atom_chunk_size, &
1099 378 : explicit=onedft_atom_chunk_size_explicit)
1100 : CALL section_vals_val_get( &
1101 : gauxc_functional_section, &
1102 : "PERIODIC_REFERENCE", &
1103 378 : l_val=periodic_reference)
1104 : CALL section_vals_val_get( &
1105 : gauxc_functional_section, &
1106 : "MOLECULAR_VIRIAL", &
1107 378 : l_val=molecular_virial)
1108 : CALL section_vals_val_get( &
1109 : gauxc_functional_section, &
1110 : "MOLECULAR_VIRIAL_DEBUG", &
1111 378 : l_val=molecular_virial_debug)
1112 : CALL section_vals_val_get( &
1113 : gauxc_functional_section, &
1114 : "MOLECULAR_VIRIAL_DEBUG_DX", &
1115 378 : r_val=molecular_virial_debug_dx)
1116 : CALL section_vals_val_get( &
1117 : gauxc_functional_section, &
1118 : "LB_EXECUTION_SPACE", &
1119 378 : c_val=lb_exec_space)
1120 : CALL section_vals_val_get( &
1121 : gauxc_functional_section, &
1122 : "INT_EXECUTION_SPACE", &
1123 378 : c_val=int_exec_space)
1124 : CALL section_vals_val_get( &
1125 : gauxc_functional_section, &
1126 : "LWD_KERNEL", &
1127 378 : c_val=lwd_kernel)
1128 : CALL section_vals_val_get( &
1129 : gauxc_functional_section, &
1130 : "SKALA_RUNTIME", &
1131 378 : c_val=skala_runtime)
1132 : CALL section_vals_val_get( &
1133 : gauxc_functional_section, &
1134 : "ONEDFT_GRADIENT_RUNTIME", &
1135 378 : c_val=onedft_gradient_runtime)
1136 : CALL section_vals_val_get( &
1137 : gauxc_functional_section, &
1138 : "OUTPUT_PATH", &
1139 378 : c_val=output_path)
1140 :
1141 378 : model_key = ADJUSTL(model_name)
1142 378 : CALL uppercase(model_key)
1143 378 : skala_runtime_key = ADJUSTL(skala_runtime)
1144 378 : CALL uppercase(skala_runtime_key)
1145 378 : onedft_gradient_runtime_key = ADJUSTL(onedft_gradient_runtime)
1146 378 : CALL uppercase(onedft_gradient_runtime_key)
1147 378 : use_onedft = (TRIM(model_key) /= "" .AND. TRIM(model_key) /= "NONE")
1148 378 : use_skala_model = (INDEX(TRIM(model_key), "SKALA") > 0)
1149 378 : IF (gapw_pseudopotentials .AND. .NOT. use_onedft) THEN
1150 : CALL cp_abort(__LOCATION__, &
1151 : "GauXC with METHOD GAPW and pseudopotentials is supported only for "// &
1152 : "OneDFT/SKALA-style models that replace the molecular XC term. "// &
1153 : "Use POTENTIAL ALL for local/semi-local GauXC GAPW validation or METHOD GPW "// &
1154 0 : "with pseudopotentials.")
1155 : END IF
1156 378 : IF (gapw_paw_pseudopotentials .AND. use_onedft) THEN
1157 : CALL cp_abort(__LOCATION__, &
1158 : "GauXC OneDFT/SKALA with METHOD GAPW and GTH/ECP pseudopotentials supports "// &
1159 : "only non-PAW regular-grid kinds, for example kinds treated through GPW_TYPE. "// &
1160 : "PAW/one-center GAPW pseudopotential kinds need a dedicated SKALA-consistent "// &
1161 0 : "one-center density design.")
1162 : END IF
1163 378 : IF (gapw_pseudopotentials .AND. use_onedft .AND. para_env%mepos == 0 .AND. &
1164 : ASSOCIATED(scf_env)) THEN
1165 2 : IF (scf_env%iter_count == 1) THEN
1166 : CALL cp_warn( &
1167 : __LOCATION__, &
1168 : "GauXC OneDFT/SKALA with METHOD GAPW and pseudopotentials evaluates the XC term "// &
1169 : "directly on the molecular AO/valence density. CP2K's GAPW one-center XC "// &
1170 2 : "correction is not used; METHOD GAPW_XC with GauXC remains unsupported.")
1171 : END IF
1172 : END IF
1173 378 : IF (device_runtime_fill_fraction <= 0.0_dp .OR. device_runtime_fill_fraction > 1.0_dp) THEN
1174 : CALL cp_abort(__LOCATION__, &
1175 0 : "GAUXC%DEVICE_RUNTIME_FILL_FRACTION must be > 0 and <= 1.")
1176 : END IF
1177 378 : IF (onedft_atom_chunk_size < -1) THEN
1178 : CALL cp_abort(__LOCATION__, &
1179 0 : "GAUXC%ONEDFT_ATOM_CHUNK_SIZE must be -1, zero, or positive.")
1180 : END IF
1181 378 : IF (molecular_virial_debug) THEN
1182 0 : IF (molecular_virial_debug_dx <= 0.0_dp) THEN
1183 : CALL cp_abort(__LOCATION__, &
1184 0 : "GauXC MOLECULAR_VIRIAL_DEBUG_DX must be positive.")
1185 : END IF
1186 0 : molecular_virial = .TRUE.
1187 : END IF
1188 378 : IF (gapw_pseudopotentials .AND. use_onedft .AND. &
1189 : (calculate_forces .OR. molecular_virial)) THEN
1190 : CALL cp_abort(__LOCATION__, &
1191 : "GauXC OneDFT/SKALA with METHOD GAPW and pseudopotentials currently "// &
1192 : "supports energies only. Nuclear gradients and molecular virials need a "// &
1193 0 : "dedicated derivative of the molecular AO/valence-density XC path.")
1194 : END IF
1195 : CALL ensure_gauxc_periodic_reference_scope( &
1196 378 : dft_control, cell, qs_kind_set, do_kpoints, periodic_reference)
1197 378 : IF (is_periodic .AND. periodic_reference .AND. para_env%mepos == 0) THEN
1198 219 : IF (ASSOCIATED(scf_env)) THEN
1199 219 : IF (scf_env%iter_count == 1) THEN
1200 : CALL cp_warn( &
1201 : __LOCATION__, &
1202 : "GAUXC%PERIODIC_REFERENCE uses GauXC molecular quadrature for isolated validation "// &
1203 30 : "cells. Compact periodic materials require a dedicated periodic GauXC interface.")
1204 : END IF
1205 : END IF
1206 : END IF
1207 378 : IF (dft_control%qs_control%gapw_xc .AND. use_onedft) THEN
1208 : CALL cp_abort(__LOCATION__, &
1209 : "GauXC OneDFT/SKALA with METHOD GAPW_XC is not implemented. "// &
1210 0 : "The GAPW_XC one-center XC correction must be evaluated by the same GauXC model.")
1211 : END IF
1212 378 : IF (use_onedft) THEN
1213 354 : IF (has_nlcc(qs_kind_set)) THEN
1214 : CALL cp_abort(__LOCATION__, &
1215 : "GauXC OneDFT/SKALA with NLCC pseudopotentials is not implemented. "// &
1216 0 : "The frozen core density would need a SKALA-consistent feature definition.")
1217 : END IF
1218 : END IF
1219 354 : IF (use_onedft) THEN
1220 : CALL set_gauxc_onedft_atom_chunk_env( &
1221 354 : onedft_atom_chunk_size, onedft_atom_chunk_size_explicit)
1222 354 : IF (.NOT. grid_explicit) grid_type = "SUPERFINE"
1223 354 : IF (.NOT. pruning_explicit) pruning_scheme = "UNPRUNED"
1224 :
1225 354 : grid_key = ADJUSTL(grid_type)
1226 354 : pruning_key = ADJUSTL(pruning_scheme)
1227 354 : CALL uppercase(grid_key)
1228 354 : CALL uppercase(pruning_key)
1229 354 : IF (use_skala_model .AND. (calculate_forces .OR. molecular_virial) .AND. &
1230 : (TRIM(grid_key) /= "SUPERFINE" .OR. TRIM(pruning_key) /= "UNPRUNED")) THEN
1231 : CALL cp_warn( &
1232 : __LOCATION__, &
1233 : "GauXC OneDFT/SKALA nuclear gradients are sensitive to the GauXC molecular grid. "// &
1234 0 : "Use GRID SUPERFINE and PRUNING_SCHEME UNPRUNED for quantitative force checks.")
1235 : END IF
1236 354 : IF (TRIM(model_key) == "SKALA") THEN
1237 26 : CALL GET_ENVIRONMENT_VARIABLE("GAUXC_SKALA_MODEL", model_name, STATUS=env_status)
1238 26 : IF (env_status /= 0 .OR. LEN_TRIM(model_name) == 0) THEN
1239 0 : CPABORT("MODEL SKALA requires the GAUXC_SKALA_MODEL environment variable")
1240 : END IF
1241 : END IF
1242 : END IF
1243 756 : SELECT CASE (TRIM(skala_runtime_key))
1244 : CASE ("AUTO")
1245 378 : use_self_runtime = use_skala_model .AND. para_env%num_pe > 1 .AND. nspins > 1
1246 : CASE ("MPI")
1247 0 : use_self_runtime = .FALSE.
1248 : CASE ("SELF")
1249 0 : use_self_runtime = use_skala_model .AND. para_env%num_pe > 1
1250 : CASE DEFAULT
1251 378 : CALL cp_abort(__LOCATION__, "Unknown GAUXC%SKALA_RUNTIME value.")
1252 : END SELECT
1253 26 : IF (.NOT. use_skala_model) use_self_runtime = .FALSE.
1254 756 : SELECT CASE (TRIM(onedft_gradient_runtime_key))
1255 : CASE ("AUTO", "SELF")
1256 378 : use_gradient_mpi_runtime = .FALSE.
1257 : use_gradient_self_runtime = calculate_forces .AND. use_onedft .AND. &
1258 378 : para_env%num_pe > 1 .AND. .NOT. use_self_runtime
1259 : CASE ("MPI")
1260 0 : use_gradient_mpi_runtime = calculate_forces .AND. use_onedft .AND. para_env%num_pe > 1
1261 0 : use_gradient_self_runtime = .FALSE.
1262 : CASE DEFAULT
1263 378 : CALL cp_abort(__LOCATION__, "Unknown GAUXC%ONEDFT_GRADIENT_RUNTIME value.")
1264 : END SELECT
1265 378 : IF (.NOT. use_onedft) THEN
1266 : use_gradient_mpi_runtime = .FALSE.
1267 : use_gradient_self_runtime = .FALSE.
1268 : END IF
1269 : IF (use_skala_model .AND. para_env%num_pe > 1 .AND. .NOT. use_self_runtime .AND. &
1270 378 : para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
1271 9 : IF (scf_env%iter_count == 1) THEN
1272 : CALL cp_warn( &
1273 : __LOCATION__, &
1274 : "GAUXC%SKALA_RUNTIME uses the MPI communicator for energy/VXC. "// &
1275 : "SKALA Torch atom chunks can be distributed across MPI ranks; "// &
1276 9 : "set GAUXC_ONEDFT_DISTRIBUTED_TORCH=0 to force rank-0 Torch inference.")
1277 : END IF
1278 : END IF
1279 :
1280 : gauxc_mol = gauxc_create_molecule( &
1281 : particle_set, &
1282 378 : gauxc_status)
1283 378 : CALL gauxc_check_status(gauxc_status)
1284 : gauxc_basis = gauxc_create_basisset( &
1285 : qs_kind_set, &
1286 : particle_set, &
1287 378 : gauxc_status)
1288 378 : CALL gauxc_check_status(gauxc_status)
1289 378 : hdf5_output = (TRIM(output_path) /= "")
1290 378 : write_hdf5_output = hdf5_output .AND. para_env%mepos == 0
1291 0 : IF (write_hdf5_output .AND. ASSOCIATED(scf_env)) THEN
1292 0 : write_hdf5_output = scf_env%iter_count == 1
1293 : END IF
1294 0 : IF (write_hdf5_output) THEN
1295 : CALL gauxc_write_molecule_hdf5( &
1296 : gauxc_mol, &
1297 : output_path, &
1298 : "molecule.h5", &
1299 : "molecule", &
1300 0 : gauxc_status)
1301 0 : CALL gauxc_check_status(gauxc_status)
1302 : CALL gauxc_write_basisset_hdf5( &
1303 : gauxc_basis, &
1304 : output_path, &
1305 : "basisset.h5", &
1306 : "basisset", &
1307 0 : gauxc_status)
1308 0 : CALL gauxc_check_status(gauxc_status)
1309 : END IF
1310 : use_fd_gradient = dft_control%qs_control%gapw .AND. calculate_forces .AND. &
1311 378 : (use_skala_model .OR. gauxc_basis%max_l > 3)
1312 : IF (use_fd_gradient) use_gradient_mpi_runtime = .FALSE.
1313 378 : IF (use_gradient_mpi_runtime .AND. use_self_runtime) THEN
1314 : CALL cp_abort( &
1315 : __LOCATION__, &
1316 : "GAUXC%ONEDFT_GRADIENT_RUNTIME MPI requires a GauXC MPI runtime. "// &
1317 0 : "Use GAUXC%SKALA_RUNTIME MPI or a closed-shell AUTO runtime.")
1318 : END IF
1319 378 : IF (use_fd_gradient .AND. para_env%mepos == 0) THEN
1320 : CALL cp_warn( &
1321 : __LOCATION__, &
1322 : "Using finite-difference GauXC XC gradients for METHOD GAPW with SKALA or "// &
1323 : "all-electron basis functions beyond f shells. The upstream analytical GauXC "// &
1324 0 : "gradient path is not yet reliable for this case.")
1325 : END IF
1326 378 : IF (use_self_runtime) THEN
1327 : ! SKALA currently needs a replicated molecular runtime for reproducible
1328 : ! open-shell densities across CP2K MPI ranks.
1329 : gauxc_grid_result = gauxc_create_grid( &
1330 : gauxc_mol, &
1331 : gauxc_basis, &
1332 : grid_type, &
1333 : radial_quadrature, &
1334 : pruning_scheme, &
1335 : lb_exec_space, &
1336 : batch_size, &
1337 : device_runtime_fill_fraction, &
1338 : gauxc_status, &
1339 8 : mpi_comm=mp_comm_self%get_handle())
1340 : ELSE
1341 : ! Use the QS force-evaluation communicator rather than GauXC's global
1342 : ! runtime. In mixed CDFT the diabatic states can be built on disjoint
1343 : ! MPI subgroups; using the global communicator would make GauXC wait
1344 : ! for ranks that are working on another state.
1345 : gauxc_grid_result = gauxc_create_grid( &
1346 : gauxc_mol, &
1347 : gauxc_basis, &
1348 : grid_type, &
1349 : radial_quadrature, &
1350 : pruning_scheme, &
1351 : lb_exec_space, &
1352 : batch_size, &
1353 : device_runtime_fill_fraction, &
1354 : gauxc_status, &
1355 370 : mpi_comm=para_env%get_handle())
1356 : END IF
1357 378 : CALL gauxc_check_status(gauxc_status)
1358 : gauxc_integrator_result = gauxc_create_integrator( &
1359 : TRIM(xc_fun_name), &
1360 : gauxc_grid_result, &
1361 : int_exec_space, &
1362 : lwd_kernel, &
1363 : nspins, &
1364 378 : gauxc_status)
1365 378 : CALL gauxc_check_status(gauxc_status)
1366 :
1367 378 : IF (use_gradient_self_runtime) THEN
1368 : ! Upstream GauXC does not yet support OneDFT/SKALA nuclear gradients
1369 : ! on an MPI runtime. Keep the energy/VXC path on the normal MPI
1370 : ! runtime and use an isolated runtime only for the replicated gradient.
1371 : gauxc_gradient_grid_result = gauxc_create_grid( &
1372 : gauxc_mol, &
1373 : gauxc_basis, &
1374 : grid_type, &
1375 : radial_quadrature, &
1376 : pruning_scheme, &
1377 : lb_exec_space, &
1378 : batch_size, &
1379 : device_runtime_fill_fraction, &
1380 : gauxc_status, &
1381 4 : mpi_comm=mp_comm_self%get_handle())
1382 4 : CALL gauxc_check_status(gauxc_status)
1383 : gauxc_gradient_integrator_result = gauxc_create_integrator( &
1384 : TRIM(xc_fun_name), &
1385 : gauxc_gradient_grid_result, &
1386 : int_exec_space, &
1387 : lwd_kernel, &
1388 : nspins, &
1389 4 : gauxc_status)
1390 4 : CALL gauxc_check_status(gauxc_status)
1391 : END IF
1392 :
1393 378 : IF (qs_env%run_rtp) THEN
1394 0 : CPABORT("GAUXC XC energy currently does not support real-time propagation")
1395 : END IF
1396 :
1397 378 : energy%exc = 0
1398 :
1399 378 : IF (ASSOCIATED(matrix_vxc)) CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1400 378 : CALL dbcsr_allocate_matrix_set(matrix_vxc, nspins)
1401 :
1402 756 : DO img = 1, nimages
1403 378 : IF (img > 1) THEN
1404 0 : CPABORT("UNIMPLEMENTED: Handling nimg>1 in k-point integration")
1405 : END IF
1406 378 : CALL dbcsr_to_dense(rho_ao(1, img), density_scalar)
1407 378 : CALL para_env%sum(density_scalar)
1408 378 : IF (nspins == 1) THEN
1409 : gauxc_xc_result = gauxc_compute_xc( &
1410 : gauxc_integrator_result, &
1411 : density_scalar, &
1412 : nspins=nspins, &
1413 : status=gauxc_status, &
1414 356 : model=TRIM(model_name))
1415 356 : CALL gauxc_check_status(gauxc_status)
1416 356 : IF (calculate_forces) THEN
1417 4 : IF (use_fd_gradient) THEN
1418 : CALL gauxc_xc_gradient_fd( &
1419 : particle_set, qs_kind_set, density_scalar, nspins, model_name, &
1420 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
1421 : lb_exec_space, int_exec_space, lwd_kernel, batch_size, &
1422 : device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
1423 0 : exc_grad%exc_grad)
1424 4 : ELSE IF (use_gradient_self_runtime) THEN
1425 : exc_grad = gauxc_compute_xc_gradient( &
1426 : gauxc_gradient_integrator_result, &
1427 : density_scalar, &
1428 : nspins=nspins, &
1429 : natom=natom, &
1430 : status=gauxc_status, &
1431 4 : model=TRIM(model_name))
1432 : ELSE
1433 : exc_grad = gauxc_compute_xc_gradient( &
1434 : gauxc_integrator_result, &
1435 : density_scalar, &
1436 : nspins=nspins, &
1437 : natom=natom, &
1438 : status=gauxc_status, &
1439 0 : model=TRIM(model_name))
1440 : END IF
1441 4 : CALL gauxc_check_status(gauxc_status)
1442 : CALL add_gauxc_gradient_to_force( &
1443 : exc_grad%exc_grad, &
1444 : force, &
1445 : atomic_kind_set, &
1446 4 : para_env)
1447 4 : IF (molecular_virial) THEN
1448 0 : CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1449 : END IF
1450 4 : IF (molecular_virial_debug) THEN
1451 : CALL debug_gauxc_molecular_virial( &
1452 : exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, nspins, &
1453 : model_name, xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
1454 : lb_exec_space, int_exec_space, lwd_kernel, batch_size, &
1455 0 : device_runtime_fill_fraction, molecular_virial_debug_dx, para_env)
1456 : END IF
1457 4 : DEALLOCATE (exc_grad%exc_grad)
1458 : END IF
1459 : ELSE
1460 22 : CPASSERT(nspins == 2)
1461 : ! In here:
1462 : ! scalar <- rho_ao(1, :) + rho_ao(2, :)
1463 : ! zeta <- rho_ao(1, :) - rho_ao(2, :)
1464 22 : CALL dbcsr_to_dense(rho_ao(2, img), density_zeta)
1465 22 : CALL para_env%sum(density_zeta)
1466 : ! Do NOT reorder the following lines!
1467 5330 : density_scalar(:, :) = density_scalar(:, :) + density_zeta(:, :)
1468 : ! Factor two because the next line is evaluated after the above line.
1469 : ! We need to subtract density_zeta once to undo the above line and
1470 : ! a second time because that is what UKS requires.
1471 : ! This style lowers memory footprint.
1472 5330 : density_zeta(:, :) = density_scalar(:, :) - 2.0_dp*density_zeta(:, :)
1473 : gauxc_xc_result = gauxc_compute_xc( &
1474 : gauxc_integrator_result, &
1475 : density_scalar, &
1476 : density_zeta, &
1477 : nspins, &
1478 : gauxc_status, &
1479 22 : model=TRIM(model_name))
1480 22 : CALL gauxc_check_status(gauxc_status)
1481 22 : IF (calculate_forces) THEN
1482 0 : IF (use_fd_gradient) THEN
1483 : CALL gauxc_xc_gradient_fd( &
1484 : particle_set, qs_kind_set, density_scalar, nspins, model_name, &
1485 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
1486 : lb_exec_space, int_exec_space, lwd_kernel, batch_size, &
1487 : device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
1488 0 : exc_grad%exc_grad, density_zeta=density_zeta)
1489 0 : ELSE IF (use_gradient_self_runtime) THEN
1490 : exc_grad = gauxc_compute_xc_gradient( &
1491 : gauxc_gradient_integrator_result, &
1492 : density_scalar, &
1493 : density_zeta, &
1494 : nspins, &
1495 : natom, &
1496 : gauxc_status, &
1497 0 : model=TRIM(model_name))
1498 : ELSE
1499 : exc_grad = gauxc_compute_xc_gradient( &
1500 : gauxc_integrator_result, &
1501 : density_scalar, &
1502 : density_zeta, &
1503 : nspins, &
1504 : natom, &
1505 : gauxc_status, &
1506 0 : model=TRIM(model_name))
1507 : END IF
1508 0 : CALL gauxc_check_status(gauxc_status)
1509 : CALL add_gauxc_gradient_to_force( &
1510 : exc_grad%exc_grad, &
1511 : force, &
1512 : atomic_kind_set, &
1513 0 : para_env)
1514 0 : IF (molecular_virial) THEN
1515 0 : CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1516 : END IF
1517 0 : IF (molecular_virial_debug) THEN
1518 : CALL debug_gauxc_molecular_virial( &
1519 : exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, nspins, &
1520 : model_name, xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
1521 : lb_exec_space, int_exec_space, lwd_kernel, batch_size, &
1522 : device_runtime_fill_fraction, molecular_virial_debug_dx, para_env, &
1523 0 : density_zeta=density_zeta)
1524 : END IF
1525 0 : DEALLOCATE (exc_grad%exc_grad)
1526 : END IF
1527 : END IF
1528 :
1529 378 : energy%exc = energy%exc + gauxc_xc_result%exc
1530 :
1531 756 : IF (nspins == 1) THEN
1532 356 : IF (img == 1) THEN
1533 356 : matrix_vxc(1) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(1, img))
1534 : ELSE
1535 0 : CPABORT("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1536 : END IF
1537 : ELSE
1538 22 : CPASSERT(nspins == 2)
1539 : ! Transform derivatives from total/spin density back to alpha/beta channels.
1540 22 : vxc_zeta_tmp = dense_to_dbcsr(gauxc_xc_result%vxc_zeta, rho_ao(1, img))
1541 22 : IF (img == 1) THEN
1542 66 : DO ispin = 1, 2
1543 44 : matrix_vxc(ispin) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(ispin, 1))
1544 : CALL dbcsr_add( &
1545 : matrix_vxc(ispin)%matrix, &
1546 : vxc_zeta_tmp%matrix, &
1547 : 1.0_dp, &
1548 : ! 1.0 for ispin==1, -1.0 for ispin==2
1549 66 : 1.0_dp - REAL(ispin - 1, dp)*2.0_dp)
1550 : END DO
1551 : ELSE
1552 0 : CPABORT("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1553 : END IF
1554 22 : CALL dbcsr_release(vxc_zeta_tmp%matrix)
1555 22 : DEALLOCATE (vxc_zeta_tmp%matrix)
1556 : END IF
1557 : END DO
1558 :
1559 378 : DEALLOCATE (density_scalar)
1560 378 : IF (ALLOCATED(density_zeta)) DEALLOCATE (density_zeta)
1561 378 : DEALLOCATE (gauxc_xc_result%vxc_scalar)
1562 378 : IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
1563 :
1564 378 : CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
1565 778 : DO ispin = 1, nspins
1566 778 : CALL dbcsr_finalize(matrix_vxc(ispin)%matrix)
1567 : END DO
1568 :
1569 378 : IF (use_gradient_self_runtime) THEN
1570 4 : CALL gauxc_destroy_integrator(gauxc_gradient_integrator_result, gauxc_status)
1571 4 : CALL gauxc_check_status(gauxc_status)
1572 4 : CALL gauxc_destroy_grid(gauxc_gradient_grid_result, gauxc_status)
1573 4 : CALL gauxc_check_status(gauxc_status)
1574 : END IF
1575 378 : CALL gauxc_destroy_integrator(gauxc_integrator_result, gauxc_status)
1576 378 : CALL gauxc_check_status(gauxc_status)
1577 378 : CALL gauxc_destroy_grid(gauxc_grid_result, gauxc_status)
1578 378 : CALL gauxc_check_status(gauxc_status)
1579 378 : CALL gauxc_destroy_basisset(gauxc_basis, gauxc_status)
1580 378 : CALL gauxc_check_status(gauxc_status)
1581 378 : CALL gauxc_destroy_molecule(gauxc_mol, gauxc_status)
1582 378 : CALL gauxc_check_status(gauxc_status)
1583 :
1584 756 : END SUBROUTINE apply_gauxc
1585 :
1586 : END MODULE xc_gauxc_functional
|
