(git:b6ef100)
Loading...
Searching...
No Matches
xc_gauxc_functional.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8#ifdef __GAUXC
9#include "gauxc/gauxc_config.f"
10#endif
11
12MODULE xc_gauxc_functional
13 USE atomic_kind_types, ONLY: atomic_kind_type,&
14 get_atomic_kind
15 USE cell_types, ONLY: cell_type
16 USE cp_control_types, ONLY: dft_control_type
17 USE cp_dbcsr_api, ONLY: dbcsr_add,&
18 dbcsr_create,&
19 dbcsr_finalize,&
20 dbcsr_get_info,&
21 dbcsr_p_type,&
22 dbcsr_release
23 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
24 dbcsr_deallocate_matrix_set
25 USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
26 USE external_potential_types, ONLY: all_potential_type,&
27 gth_potential_type,&
28 sgp_potential_type
29 USE input_constants, ONLY: xc_vdw_fun_nonloc
30 USE input_section_types, ONLY: section_vals_get_subs_vals,&
31 section_vals_get_subs_vals2,&
32 section_vals_type,&
33 section_vals_val_get
34 USE iso_c_binding, ONLY: c_char,&
35 c_double,&
36 c_int,&
37 c_null_char
38 USE kinds, ONLY: default_path_length,&
39 default_string_length,&
40 dp
41 USE message_passing, ONLY: mp_comm_self,&
42 mp_para_env_type
43 USE particle_types, ONLY: particle_type
44 USE qs_energy_types, ONLY: qs_energy_type
45 USE qs_environment_types, ONLY: get_qs_env,&
46 qs_environment_type
47 USE qs_force_types, ONLY: qs_force_type
48 USE qs_kind_types, ONLY: get_qs_kind,&
49 has_nlcc,&
50 qs_kind_type
51 USE qs_ks_types, ONLY: qs_ks_env_type,&
52 set_ks_env
53 USE qs_rho_types, ONLY: qs_rho_get,&
54 qs_rho_type
55 USE qs_scf_types, ONLY: qs_scf_env_type
56 USE string_utilities, ONLY: uppercase
57 USE xc_gauxc_cache, ONLY: cp_gauxc_cache_params,&
58 cp_gauxc_cache_type,&
59 gauxc_cache_init
60 USE xc_gauxc_interface, ONLY: &
61 cp_gauxc_basisset_type, cp_gauxc_grid_type, cp_gauxc_integrator_type, &
62 cp_gauxc_molecule_type, cp_gauxc_status_type, cp_gauxc_xc_gradient_type, cp_gauxc_xc_type, &
63 gauxc_check_status, gauxc_compute_xc, gauxc_compute_xc_gradient, gauxc_create_basisset, &
64 gauxc_create_grid, gauxc_create_integrator, gauxc_create_molecule, gauxc_destroy_basisset, &
65 gauxc_destroy_grid, gauxc_destroy_integrator, gauxc_destroy_molecule, &
66 gauxc_write_basisset_hdf5, gauxc_write_molecule_hdf5
67 USE xc_input_constants, ONLY: skala_gapw_paw_one_center
68 USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
69#include "../base/base_uses.f90"
70
71 IMPLICIT NONE
72
73 PRIVATE
74
75 LOGICAL, PARAMETER :: debug_this_module = .true.
76 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_gauxc_functional'
77
78 PUBLIC :: apply_gauxc, gauxc_gapw_all_all_electron, gauxc_gapw_all_pseudopotentials, &
79 gauxc_gapw_has_paw_pseudopotentials, gauxc_gapw_has_pseudopotentials, skala_info, &
80 xc_section_uses_gauxc
81
82 INTERFACE
83 INTEGER(c_int) FUNCTION c_setenv(name, value, overwrite) BIND(C, name="setenv")
84 IMPORT :: c_char, c_int
85 CHARACTER(KIND=c_char), DIMENSION(*), INTENT(IN) :: name, value
86 INTEGER(c_int), VALUE :: overwrite
87 END FUNCTION c_setenv
88
89 INTEGER(c_int) FUNCTION c_unsetenv(name) BIND(C, name="unsetenv")
90 IMPORT :: c_char, c_int
91 CHARACTER(KIND=c_char), DIMENSION(*), INTENT(IN) :: name
92 END FUNCTION c_unsetenv
93 END INTERFACE
94
95CONTAINS
96
97! **************************************************************************************************
98!> \brief Return whether every GAPW kind uses an all-electron potential.
99!> \param qs_kind_set ...
100!> \return ...
101! **************************************************************************************************
102 FUNCTION gauxc_gapw_all_all_electron(qs_kind_set) RESULT(all_all_electron)
103 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
104 LOGICAL :: all_all_electron
105
106 INTEGER :: ikind
107 TYPE(all_potential_type), POINTER :: all_potential
108
109 cpassert(ASSOCIATED(qs_kind_set))
110
111 all_all_electron = .true.
112 DO ikind = 1, SIZE(qs_kind_set)
113 NULLIFY (all_potential)
114 CALL get_qs_kind(qs_kind_set(ikind), all_potential=all_potential)
115 IF (.NOT. ASSOCIATED(all_potential)) THEN
116 all_all_electron = .false.
117 EXIT
118 END IF
119 END DO
120
121 END FUNCTION gauxc_gapw_all_all_electron
122
123! **************************************************************************************************
124!> \brief Set the GauXC Skala atom chunk environment knob when the CP2K keyword is explicit.
125!> \param atom_chunk_size ...
126!> \param is_explicit ...
127! **************************************************************************************************
128 SUBROUTINE set_gauxc_model_atom_chunk_env(atom_chunk_size, is_explicit)
129 INTEGER, INTENT(IN) :: atom_chunk_size
130 LOGICAL, INTENT(IN) :: is_explicit
131
132 CHARACTER(LEN=32) :: chunk_value
133 INTEGER(c_int) :: ierr
134
135 IF (.NOT. is_explicit) RETURN
136
137 IF (atom_chunk_size < 0) THEN
138 ierr = c_unsetenv("GAUXC_ONEDFT_ATOM_CHUNK_SIZE"//c_null_char)
139 ELSE
140 WRITE (chunk_value, '(I0)') atom_chunk_size
141 ierr = c_setenv( &
142 "GAUXC_ONEDFT_ATOM_CHUNK_SIZE"//c_null_char, &
143 trim(chunk_value)//c_null_char, &
144 1_c_int)
145 END IF
146 IF (ierr /= 0_c_int) THEN
147 CALL cp_abort(__location__, &
148 "Could not set GAUXC_ONEDFT_ATOM_CHUNK_SIZE for GauXC Skala.")
149 END IF
150 END SUBROUTINE set_gauxc_model_atom_chunk_env
151
152! **************************************************************************************************
153!> \brief ...
154!> \param dbcsr_mat ...
155!> \param dense_mat ...
156!> \param para_env ...
157! **************************************************************************************************
158 SUBROUTINE dbcsr_to_dense(dbcsr_mat, dense_mat, para_env)
159 USE cp_dbcsr_api, ONLY: dbcsr_distribution_get, dbcsr_distribution_type, dbcsr_get_info, &
160 dbcsr_get_matrix_type, dbcsr_get_readonly_block_p, &
161 dbcsr_get_stored_coordinates, dbcsr_type_antisymmetric, &
162 dbcsr_type_symmetric
163 TYPE(dbcsr_p_type), INTENT(IN) :: dbcsr_mat
164 REAL(c_double), ALLOCATABLE, DIMENSION(:, :), &
165 INTENT(INOUT) :: dense_mat
166 TYPE(mp_para_env_type), INTENT(IN), POINTER :: para_env
167
168 CHARACTER :: matrix_type
169 INTEGER :: col, col_end, col_start, icol, irow, mynode, nblkcols_total, nblkrows_total, &
170 ncols, nrows, numnodes, owner, row, row_end, row_start
171 INTEGER, ALLOCATABLE, DIMENSION(:) :: c_offset, r_offset
172 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
173 LOGICAL :: found
174 REAL(c_double), POINTER :: block(:, :)
175 TYPE(dbcsr_distribution_type) :: dist
176
177 CALL dbcsr_get_info(dbcsr_mat%matrix, &
178 row_blk_size=row_blk_size, &
179 col_blk_size=col_blk_size, &
180 nblkrows_total=nblkrows_total, &
181 nblkcols_total=nblkcols_total, &
182 nfullrows_total=nrows, &
183 nfullcols_total=ncols, &
184 distribution=dist)
185 CALL dbcsr_distribution_get(dist, mynode=mynode, numnodes=numnodes)
186 matrix_type = dbcsr_get_matrix_type(dbcsr_mat%matrix)
187
188 IF (.NOT. ALLOCATED(dense_mat)) THEN
189 ALLOCATE (dense_mat(nrows, ncols))
190 ELSE IF (.NOT. all(shape(dense_mat) == [nrows, ncols])) THEN
191 DEALLOCATE (dense_mat)
192 ALLOCATE (dense_mat(nrows, ncols))
193 ELSE
194 cpassert(all(shape(dense_mat) == [nrows, ncols]))
195 END IF
196 dense_mat = 0._dp
197
198 ALLOCATE (r_offset(nblkrows_total), c_offset(nblkcols_total))
199
200 r_offset(1) = 1
201 DO row = 2, nblkrows_total
202 r_offset(row) = r_offset(row - 1) + row_blk_size(row - 1)
203 END DO
204 c_offset(1) = 1
205 DO col = 2, nblkcols_total
206 c_offset(col) = c_offset(col - 1) + col_blk_size(col - 1)
207 END DO
208
209 ! Replicated DBCSR blocks must enter the following MPI sum exactly once.
210 DO irow = 1, nblkrows_total
211 DO icol = 1, nblkcols_total
212 IF (numnodes == 1 .AND. para_env%num_pe > 1 .AND. para_env%mepos /= 0) cycle
213 CALL dbcsr_get_stored_coordinates(dbcsr_mat%matrix, irow, icol, owner)
214 IF (owner /= mynode) cycle
215 CALL dbcsr_get_readonly_block_p(matrix=dbcsr_mat%matrix, row=irow, col=icol, &
216 block=block, found=found)
217 IF (.NOT. found) cycle
218 row_start = r_offset(irow)
219 row_end = row_start + row_blk_size(irow) - 1
220 col_start = c_offset(icol)
221 col_end = col_start + col_blk_size(icol) - 1
222 dense_mat(row_start:row_end, col_start:col_end) = block
223 IF (irow /= icol) THEN
224 IF (matrix_type == dbcsr_type_symmetric) THEN
225 dense_mat(col_start:col_end, row_start:row_end) = transpose(block)
226 ELSE IF (matrix_type == dbcsr_type_antisymmetric) THEN
227 dense_mat(col_start:col_end, row_start:row_end) = -transpose(block)
228 END IF
229 END IF
230 END DO
231 END DO
232
233 DEALLOCATE (r_offset, c_offset)
234
235 END SUBROUTINE dbcsr_to_dense
236
237! ******, ***********************************************************************************
238!> \brief Convert a dense symmetric matrix to a DBCSR matrix with full upper block structure.
239!> This creates all upper-triangular blocks, not just those present in a template.
240!> This is needed because GauXC computes VXC for the full dense density matrix.
241!> \param dense_mat Input dense matrix
242!> \param template_dbcsr Template DBCSR matrix for distribution and block sizes
243!> \return dbcsr_mat Output DBCSR matrix with full upper block structure
244! **************************************************************************************************
245 FUNCTION dense_to_dbcsr(dense_mat, template_dbcsr) RESULT(dbcsr_mat)
246 USE cp_dbcsr_api, ONLY: &
247 dbcsr_create, &
248 dbcsr_distribution_get, &
249 dbcsr_distribution_type, &
250 dbcsr_finalize, &
251 dbcsr_get_info, &
252 dbcsr_get_stored_coordinates, &
253 dbcsr_init_p, &
254 dbcsr_put_block, &
255 dbcsr_release, &
256 dbcsr_type_symmetric, &
257 dbcsr_work_create
258 REAL(c_double), DIMENSION(:, :), INTENT(IN) :: dense_mat
259 TYPE(dbcsr_p_type), INTENT(IN) :: template_dbcsr
260 TYPE(dbcsr_p_type) :: dbcsr_mat
261
262 INTEGER :: col, icol, irow, mynode, nblkcols_total, &
263 nblkrows_total, ncols, nrows, owner, &
264 row
265 INTEGER, ALLOCATABLE, DIMENSION(:) :: c_offset, r_offset
266 INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
267 TYPE(dbcsr_distribution_type) :: dist
268
269 CALL dbcsr_get_info(template_dbcsr%matrix, &
270 row_blk_size=row_blk_size, &
271 col_blk_size=col_blk_size, &
272 nblkrows_total=nblkrows_total, &
273 nblkcols_total=nblkcols_total, &
274 nfullrows_total=nrows, &
275 nfullcols_total=ncols, &
276 distribution=dist)
277 CALL dbcsr_distribution_get(dist, mynode=mynode)
278
279 cpassert(nrows == SIZE(dense_mat, 1))
280 cpassert(ncols == SIZE(dense_mat, 2))
281
282 CALL dbcsr_init_p(dbcsr_mat%matrix)
283 CALL dbcsr_create(dbcsr_mat%matrix, &
284 template=template_dbcsr%matrix, &
285 name="VXC from GauXC (dense)", &
286 matrix_type=dbcsr_type_symmetric)
287 CALL dbcsr_work_create(dbcsr_mat%matrix, work_mutable=.true.)
288
289 ALLOCATE (r_offset(nblkrows_total), c_offset(nblkcols_total))
290
291 r_offset(1) = 1
292 DO row = 2, nblkrows_total
293 r_offset(row) = r_offset(row - 1) + row_blk_size(row - 1)
294 END DO
295 c_offset(1) = 1
296 DO col = 2, nblkcols_total
297 c_offset(col) = c_offset(col - 1) + col_blk_size(col - 1)
298 END DO
299
300 DO irow = 1, nblkrows_total
301 DO icol = 1, nblkcols_total
302 IF (irow > icol) cycle
303 CALL dbcsr_get_stored_coordinates(dbcsr_mat%matrix, irow, icol, owner)
304 IF (owner /= mynode) cycle
305 CALL dbcsr_put_block(dbcsr_mat%matrix, irow, icol, &
306 0.5_dp*( &
307 dense_mat(r_offset(irow):r_offset(irow) + row_blk_size(irow) - 1, &
308 c_offset(icol):c_offset(icol) + col_blk_size(icol) - 1) + &
309 transpose(dense_mat(r_offset(icol):r_offset(icol) + row_blk_size(icol) - 1, &
310 c_offset(irow):c_offset(irow) + col_blk_size(irow) - 1))))
311 END DO
312 END DO
313
314 CALL dbcsr_finalize(dbcsr_mat%matrix)
315
316 DEALLOCATE (r_offset, c_offset)
317
318 END FUNCTION dense_to_dbcsr
319
320! **************************************************************************************************
321!> \brief ...
322!> \param xc_section ...
323!> \return ...
324! **************************************************************************************************
325 FUNCTION get_gauxc_functional(xc_section) RESULT(gauxc_functional_section)
326 TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
327 TYPE(section_vals_type), POINTER :: gauxc_functional_section
328
329 INTEGER :: ifun
330 TYPE(section_vals_type), POINTER :: functionals, xc_fun
331
332 NULLIFY (gauxc_functional_section)
333
334 functionals => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
335 IF (.NOT. ASSOCIATED(functionals)) THEN
336 cpabort("XC_FUNCTIONAL section not found")
337 END IF
338
339 ifun = 0
340 DO
341 ifun = ifun + 1
342 xc_fun => section_vals_get_subs_vals2(functionals, i_section=ifun)
343 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
344 IF (xc_fun%section%name /= "GAUXC" .OR. ifun > 1) THEN
345 cpabort("GauXC functionals are mutually exclusive with any other functional.")
346 END IF
347 gauxc_functional_section => xc_fun
348 END DO
349
350 IF (.NOT. ASSOCIATED(gauxc_functional_section)) THEN
351 cpabort("No XC functional found in XC_FUNCTIONAL section")
352 END IF
353 END FUNCTION get_gauxc_functional
354
355! **************************************************************************************************
356!> \brief ...
357!> \param xc_section ...
358!> \return ...
359! **************************************************************************************************
360 FUNCTION xc_section_uses_gauxc(xc_section) RESULT(uses_gauxc)
361 TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
362 LOGICAL :: uses_gauxc
363
364 INTEGER :: ifun
365 TYPE(section_vals_type), POINTER :: functionals, xc_fun
366
367 uses_gauxc = .false.
368 IF (.NOT. ASSOCIATED(xc_section)) RETURN
369
370 functionals => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
371 IF (.NOT. ASSOCIATED(functionals)) RETURN
372
373 ifun = 0
374 DO
375 ifun = ifun + 1
376 xc_fun => section_vals_get_subs_vals2(functionals, i_section=ifun)
377 IF (.NOT. ASSOCIATED(xc_fun)) EXIT
378 IF (xc_fun%section%name == "GAUXC") THEN
379 uses_gauxc = .true.
380 EXIT
381 END IF
382 END DO
383
384 END FUNCTION xc_section_uses_gauxc
385
386! **************************************************************************************************
387!> \brief Return whether GauXC GAPW mode sees pseudopotential kinds.
388!> \param qs_kind_set ...
389!> \return ...
390! **************************************************************************************************
391 FUNCTION gauxc_gapw_has_pseudopotentials(qs_kind_set) RESULT(has_pseudopotentials)
392 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
393 LOGICAL :: has_pseudopotentials
394
395 INTEGER :: ikind
396 TYPE(gth_potential_type), POINTER :: gth_potential
397 TYPE(sgp_potential_type), POINTER :: sgp_potential
398
399 cpassert(ASSOCIATED(qs_kind_set))
400
401 has_pseudopotentials = .false.
402 DO ikind = 1, SIZE(qs_kind_set)
403 NULLIFY (gth_potential, sgp_potential)
404 CALL get_qs_kind(qs_kind_set(ikind), &
405 gth_potential=gth_potential, &
406 sgp_potential=sgp_potential)
407 IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
408 has_pseudopotentials = .true.
409 EXIT
410 END IF
411 END DO
412
413 END FUNCTION gauxc_gapw_has_pseudopotentials
414
415! **************************************************************************************************
416!> \brief Return whether every GAPW kind uses a pseudopotential.
417!> \param qs_kind_set ...
418!> \return ...
419! **************************************************************************************************
420 FUNCTION gauxc_gapw_all_pseudopotentials(qs_kind_set) RESULT(all_pseudopotentials)
421 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
422 LOGICAL :: all_pseudopotentials
423
424 INTEGER :: ikind
425 TYPE(gth_potential_type), POINTER :: gth_potential
426 TYPE(sgp_potential_type), POINTER :: sgp_potential
427
428 cpassert(ASSOCIATED(qs_kind_set))
429
430 all_pseudopotentials = .true.
431 DO ikind = 1, SIZE(qs_kind_set)
432 NULLIFY (gth_potential, sgp_potential)
433 CALL get_qs_kind(qs_kind_set(ikind), &
434 gth_potential=gth_potential, &
435 sgp_potential=sgp_potential)
436 IF (.NOT. (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential))) THEN
437 all_pseudopotentials = .false.
438 EXIT
439 END IF
440 END DO
441
442 END FUNCTION gauxc_gapw_all_pseudopotentials
443
444! **************************************************************************************************
445!> \brief Return whether GauXC GAPW mode sees pseudopotential one-center GAPW kinds.
446!> \param qs_kind_set ...
447!> \return ...
448! **************************************************************************************************
449 FUNCTION gauxc_gapw_has_paw_pseudopotentials(qs_kind_set) RESULT(has_paw_pseudopotentials)
450 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
451 LOGICAL :: has_paw_pseudopotentials
452
453 INTEGER :: ikind
454 LOGICAL :: paw_atom
455 TYPE(gth_potential_type), POINTER :: gth_potential
456 TYPE(sgp_potential_type), POINTER :: sgp_potential
457
458 cpassert(ASSOCIATED(qs_kind_set))
459
460 has_paw_pseudopotentials = .false.
461 DO ikind = 1, SIZE(qs_kind_set)
462 NULLIFY (gth_potential, sgp_potential)
463 CALL get_qs_kind(qs_kind_set(ikind), &
464 gth_potential=gth_potential, &
465 paw_atom=paw_atom, &
466 sgp_potential=sgp_potential)
467 IF ((ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) .AND. &
468 paw_atom) THEN
469 has_paw_pseudopotentials = .true.
470 EXIT
471 END IF
472 END DO
473
474 END FUNCTION gauxc_gapw_has_paw_pseudopotentials
475
476! **************************************************************************************************
477!> \brief Check the current periodic scope of the CP2K-GauXC bridge
478!> \param dft_control ...
479!> \param cell ...
480!> \param qs_kind_set ...
481!> \param do_kpoints ...
482!> \param periodic_reference ...
483!> \note This path keeps isolated validation cells usable under PERIODIC XYZ.
484!> It intentionally does not implement compact periodic GauXC quadrature.
485! **************************************************************************************************
486 SUBROUTINE ensure_gauxc_periodic_reference_scope( &
487 dft_control, cell, qs_kind_set, do_kpoints, periodic_reference)
488 TYPE(dft_control_type), POINTER :: dft_control
489 TYPE(cell_type), POINTER :: cell
490 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
491 LOGICAL, INTENT(IN) :: do_kpoints, periodic_reference
492
493 INTEGER :: ikind
494 LOGICAL :: is_periodic
495 TYPE(gth_potential_type), POINTER :: gth_potential
496 TYPE(sgp_potential_type), POINTER :: sgp_potential
497
498 cpassert(ASSOCIATED(dft_control))
499 cpassert(ASSOCIATED(qs_kind_set))
500
501 is_periodic = .false.
502 IF (ASSOCIATED(cell)) is_periodic = any(cell%perd /= 0)
503
504 IF (do_kpoints) THEN
505 CALL cp_abort(__location__, &
506 "GauXC currently supports only Gamma-only density matrices in CP2K. "// &
507 "Periodic k-point density matrices require a dedicated GauXC periodic interface.")
508 END IF
509 IF (dft_control%nimages /= 1) THEN
510 CALL cp_abort(__location__, &
511 "GauXC currently supports only a single AO image in CP2K. "// &
512 "Periodic neighbour-cell AO blocks require a dedicated GauXC periodic interface.")
513 END IF
514 IF (.NOT. is_periodic) RETURN
515
516 IF (.NOT. periodic_reference) THEN
517 CALL cp_abort(__location__, &
518 "Periodic GauXC calculations in CP2K require GAUXC%PERIODIC_REFERENCE T. "// &
519 "This opt-in documents that the current path is only an isolated-cell, "// &
520 "Gamma-only, single-image METHOD GPW reference path using GauXC molecular "// &
521 "quadrature, not a dedicated periodic GauXC interface.")
522 END IF
523
524 IF (.NOT. all(cell%perd == 1)) THEN
525 CALL cp_abort(__location__, &
526 "The current GauXC isolated-cell reference path supports only PERIODIC XYZ. "// &
527 "Partial periodicity requires a dedicated GauXC periodic interface.")
528 END IF
529 IF (.NOT. dft_control%qs_control%gpw) THEN
530 CALL cp_abort(__location__, &
531 "The current GauXC isolated-cell reference path is limited to METHOD GPW with GTH "// &
532 "pseudopotentials. GAPW, GAPW_XC, and other QS methods are not supported here.")
533 END IF
534
535 DO ikind = 1, SIZE(qs_kind_set)
536 NULLIFY (gth_potential, sgp_potential)
537 CALL get_qs_kind(qs_kind_set(ikind), &
538 gth_potential=gth_potential, &
539 sgp_potential=sgp_potential)
540 IF (.NOT. ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
541 CALL cp_abort(__location__, &
542 "The current GauXC isolated-cell reference path is limited to GTH pseudopotentials. "// &
543 "Use non-periodic all-electron GAPW validation for molecular GAPW cases.")
544 END IF
545 END DO
546
547 END SUBROUTINE ensure_gauxc_periodic_reference_scope
548
549! **************************************************************************************************
550!> \brief adds a replicated GauXC energy gradient to the local CP2K force accumulator
551!> \param exc_grad ...
552!> \param force ...
553!> \param atomic_kind_set ...
554!> \param para_env ...
555! **************************************************************************************************
556 SUBROUTINE add_gauxc_gradient_to_force(exc_grad, force, atomic_kind_set, para_env)
557 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: exc_grad
558 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
559 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
560 TYPE(mp_para_env_type), POINTER :: para_env
561
562 INTEGER :: ia, iatom, ikind, natom_kind
563 TYPE(atomic_kind_type), POINTER :: atomic_kind
564
565 cpassert(ASSOCIATED(force))
566 cpassert(ASSOCIATED(atomic_kind_set))
567
568 IF (para_env%mepos /= 0) RETURN
569
570 DO ikind = 1, SIZE(atomic_kind_set, 1)
571 atomic_kind => atomic_kind_set(ikind)
572 CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
573 DO ia = 1, natom_kind
574 iatom = atomic_kind%atom_list(ia)
575 force(ikind)%rho_elec(:, ia) = force(ikind)%rho_elec(:, ia) + &
576 exc_grad(3*iatom - 2:3*iatom)
577 END DO
578 END DO
579
580 END SUBROUTINE add_gauxc_gradient_to_force
581
582! **************************************************************************************************
583!> \brief compute a GauXC XC energy for diagnostic finite differences
584!> \param particle_set_eval ...
585!> \param qs_kind_set ...
586!> \param density_scalar ...
587!> \param nspins ...
588!> \param model_name ...
589!> \param xc_fun_name ...
590!> \param grid_type ...
591!> \param radial_quadrature ...
592!> \param pruning_scheme ...
593!> \param lb_exec_space ...
594!> \param int_exec_space ...
595!> \param lwd_kernel ...
596!> \param batch_size ...
597!> \param device_runtime_fill_fraction ...
598!> \param exc ...
599!> \param density_zeta ...
600! **************************************************************************************************
601 SUBROUTINE gauxc_xc_energy_for_particles( &
602 particle_set_eval, qs_kind_set, density_scalar, nspins, model_name, &
603 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
604 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, exc, density_zeta)
605 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_eval
606 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
607 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
608 INTEGER, INTENT(IN) :: nspins
609 CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
610 pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
611 INTEGER, INTENT(IN) :: batch_size
612 REAL(kind=dp), INTENT(IN) :: device_runtime_fill_fraction
613 REAL(kind=dp), INTENT(OUT) :: exc
614 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
615 OPTIONAL :: density_zeta
616
617 TYPE(cp_gauxc_basisset_type) :: gauxc_basis_fd
618 TYPE(cp_gauxc_grid_type) :: gauxc_grid_fd
619 TYPE(cp_gauxc_integrator_type) :: gauxc_integrator_fd
620 TYPE(cp_gauxc_molecule_type) :: gauxc_mol_fd
621 TYPE(cp_gauxc_status_type) :: gauxc_status
622 TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
623
624 gauxc_mol_fd = gauxc_create_molecule(particle_set_eval, gauxc_status)
625 CALL gauxc_check_status(gauxc_status)
626 gauxc_basis_fd = gauxc_create_basisset(qs_kind_set, particle_set_eval, gauxc_status)
627 CALL gauxc_check_status(gauxc_status)
628 gauxc_grid_fd = gauxc_create_grid( &
629 gauxc_mol_fd, &
630 gauxc_basis_fd, &
631 grid_type, &
632 radial_quadrature, &
633 pruning_scheme, &
634 lb_exec_space, &
635 batch_size, &
636 device_runtime_fill_fraction, &
637 gauxc_status, &
638 mpi_comm=mp_comm_self%get_handle(), &
639 force_new_runtime=.true.)
640 CALL gauxc_check_status(gauxc_status)
641 gauxc_integrator_fd = gauxc_create_integrator( &
642 trim(xc_fun_name), &
643 gauxc_grid_fd, &
644 int_exec_space, &
645 lwd_kernel, &
646 nspins, &
647 gauxc_status)
648 CALL gauxc_check_status(gauxc_status)
649
650 IF (nspins == 1) THEN
651 gauxc_xc_result = gauxc_compute_xc( &
652 gauxc_integrator_fd, &
653 density_scalar, &
654 nspins=nspins, &
655 status=gauxc_status, &
656 model=trim(model_name))
657 ELSE
658 cpassert(nspins == 2)
659 cpassert(PRESENT(density_zeta))
660 gauxc_xc_result = gauxc_compute_xc( &
661 gauxc_integrator_fd, &
662 density_scalar, &
663 density_zeta, &
664 nspins, &
665 gauxc_status, &
666 model=trim(model_name))
667 END IF
668 CALL gauxc_check_status(gauxc_status)
669 exc = gauxc_xc_result%exc
670
671 IF (ALLOCATED(gauxc_xc_result%vxc_scalar)) DEALLOCATE (gauxc_xc_result%vxc_scalar)
672 IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
673
674 CALL gauxc_destroy_integrator(gauxc_integrator_fd, gauxc_status)
675 CALL gauxc_check_status(gauxc_status)
676 CALL gauxc_destroy_grid(gauxc_grid_fd, gauxc_status)
677 CALL gauxc_check_status(gauxc_status)
678 CALL gauxc_destroy_basisset(gauxc_basis_fd, gauxc_status)
679 CALL gauxc_check_status(gauxc_status)
680 CALL gauxc_destroy_molecule(gauxc_mol_fd, gauxc_status)
681 CALL gauxc_check_status(gauxc_status)
682
683 END SUBROUTINE gauxc_xc_energy_for_particles
684
685! **************************************************************************************************
686!> \brief compute a finite-difference GauXC XC nuclear gradient at fixed density
687!> \param particle_set ...
688!> \param qs_kind_set ...
689!> \param density_scalar ...
690!> \param nspins ...
691!> \param model_name ...
692!> \param xc_fun_name ...
693!> \param grid_type ...
694!> \param radial_quadrature ...
695!> \param pruning_scheme ...
696!> \param lb_exec_space ...
697!> \param int_exec_space ...
698!> \param lwd_kernel ...
699!> \param batch_size ...
700!> \param device_runtime_fill_fraction ...
701!> \param dx ...
702!> \param para_env ...
703!> \param exc_grad ...
704!> \param density_zeta ...
705! **************************************************************************************************
706 SUBROUTINE gauxc_xc_gradient_fd( &
707 particle_set, qs_kind_set, density_scalar, nspins, model_name, &
708 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
709 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, dx, para_env, exc_grad, &
710 density_zeta)
711 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
712 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
713 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
714 INTEGER, INTENT(IN) :: nspins
715 CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
716 pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
717 INTEGER, INTENT(IN) :: batch_size
718 REAL(kind=dp), INTENT(IN) :: device_runtime_fill_fraction, dx
719 TYPE(mp_para_env_type), POINTER :: para_env
720 REAL(kind=dp), ALLOCATABLE, DIMENSION(:), &
721 INTENT(OUT) :: exc_grad
722 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
723 OPTIONAL :: density_zeta
724
725 INTEGER :: iatom, idir
726 REAL(kind=dp) :: xc_minus, xc_plus
727 TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus
728
729 cpassert(ASSOCIATED(particle_set))
730 cpassert(dx > 0.0_dp)
731
732 ALLOCATE (exc_grad(3*SIZE(particle_set)))
733 exc_grad = 0.0_dp
734
735 IF (para_env%mepos == 0) THEN
736 ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set)))
737
738 DO iatom = 1, SIZE(particle_set)
739 DO idir = 1, 3
740 particle_set_minus = particle_set
741 particle_set_plus = particle_set
742 particle_set_minus(iatom)%r(idir) = particle_set_minus(iatom)%r(idir) - dx
743 particle_set_plus(iatom)%r(idir) = particle_set_plus(iatom)%r(idir) + dx
744 IF (PRESENT(density_zeta)) THEN
745 CALL gauxc_xc_energy_for_particles( &
746 particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
747 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
748 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus, &
749 density_zeta=density_zeta)
750 CALL gauxc_xc_energy_for_particles( &
751 particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
752 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
753 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus, &
754 density_zeta=density_zeta)
755 ELSE
756 CALL gauxc_xc_energy_for_particles( &
757 particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
758 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
759 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus)
760 CALL gauxc_xc_energy_for_particles( &
761 particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
762 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
763 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus)
764 END IF
765 exc_grad(3*iatom - 3 + idir) = (xc_plus - xc_minus)/(2.0_dp*dx)
766 END DO
767 END DO
768
769 DEALLOCATE (particle_set_minus, particle_set_plus)
770 END IF
771
772 CALL para_env%bcast(exc_grad, 0)
773
774 END SUBROUTINE gauxc_xc_gradient_fd
775
776! **************************************************************************************************
777!> \brief finite-difference check of the molecular GauXC XC virial diagnostic
778!> \param exc_grad ...
779!> \param particle_set ...
780!> \param qs_kind_set ...
781!> \param density_scalar ...
782!> \param nspins ...
783!> \param model_name ...
784!> \param xc_fun_name ...
785!> \param grid_type ...
786!> \param radial_quadrature ...
787!> \param pruning_scheme ...
788!> \param lb_exec_space ...
789!> \param int_exec_space ...
790!> \param lwd_kernel ...
791!> \param batch_size ...
792!> \param device_runtime_fill_fraction ...
793!> \param dx ...
794!> \param para_env ...
795!> \param density_zeta ...
796! **************************************************************************************************
797 SUBROUTINE debug_gauxc_molecular_virial( &
798 exc_grad, particle_set, qs_kind_set, density_scalar, nspins, model_name, &
799 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
800 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, dx, para_env, density_zeta)
801 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: exc_grad
802 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
803 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
804 REAL(kind=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
805 INTEGER, INTENT(IN) :: nspins
806 CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, radial_quadrature, &
807 pruning_scheme, lb_exec_space, int_exec_space, lwd_kernel
808 INTEGER, INTENT(IN) :: batch_size
809 REAL(kind=dp), INTENT(IN) :: device_runtime_fill_fraction, dx
810 TYPE(mp_para_env_type), POINTER :: para_env
811 REAL(kind=dp), DIMENSION(:, :), INTENT(IN), &
812 OPTIONAL :: density_zeta
813
814 INTEGER :: iatom, iw
815 REAL(kind=dp) :: analytic_trace, diff_trace, &
816 numerical_trace, xc_minus, xc_plus
817 REAL(kind=dp), DIMENSION(3) :: center, displacement, grad
818 TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus
819
820 cpassert(ASSOCIATED(particle_set))
821 cpassert(SIZE(exc_grad) == 3*SIZE(particle_set))
822
823 IF (para_env%mepos /= 0) RETURN
824
825 center = 0.0_dp
826 DO iatom = 1, SIZE(particle_set)
827 center = center + particle_set(iatom)%r
828 END DO
829 center = center/real(SIZE(particle_set), dp)
830
831 ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set)))
832 particle_set_minus = particle_set
833 particle_set_plus = particle_set
834
835 analytic_trace = 0.0_dp
836 DO iatom = 1, SIZE(particle_set)
837 grad = exc_grad(3*iatom - 2:3*iatom)
838 displacement = particle_set(iatom)%r - center
839 analytic_trace = analytic_trace + dot_product(grad, displacement)
840 particle_set_minus(iatom)%r = center + (1.0_dp - dx)*displacement
841 particle_set_plus(iatom)%r = center + (1.0_dp + dx)*displacement
842 END DO
843 analytic_trace = analytic_trace/3.0_dp
844
845 IF (PRESENT(density_zeta)) THEN
846 CALL gauxc_xc_energy_for_particles( &
847 particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
848 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
849 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus, &
850 density_zeta=density_zeta)
851 CALL gauxc_xc_energy_for_particles( &
852 particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
853 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
854 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus, &
855 density_zeta=density_zeta)
856 ELSE
857 CALL gauxc_xc_energy_for_particles( &
858 particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
859 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
860 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_plus)
861 CALL gauxc_xc_energy_for_particles( &
862 particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
863 xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
864 int_exec_space, lwd_kernel, batch_size, device_runtime_fill_fraction, xc_minus)
865 END IF
866
867 numerical_trace = (xc_plus - xc_minus)/(2.0_dp*dx)/3.0_dp
868 diff_trace = analytic_trace - numerical_trace
869
870 iw = cp_logger_get_default_io_unit()
871 IF (iw > 0) THEN
872 WRITE (unit=iw, fmt="(/,T2,A,1X,ES11.4)") &
873 "GAUXC| Molecular XC virial finite-difference dx", dx
874 WRITE (unit=iw, fmt="(T2,A,3(1X,ES19.11))") &
875 "GAUXC| Molecular XC virial FD 1/3 Trace", &
876 analytic_trace, numerical_trace, diff_trace
877 END IF
878
879 DEALLOCATE (particle_set_minus, particle_set_plus)
880
881 END SUBROUTINE debug_gauxc_molecular_virial
882
883! **************************************************************************************************
884!> \brief prints a force-based molecular XC virial diagnostic from GauXC gradients
885!> \param exc_grad ...
886!> \param particle_set ...
887!> \param para_env ...
888! **************************************************************************************************
889 SUBROUTINE print_gauxc_molecular_virial(exc_grad, particle_set, para_env)
890 REAL(kind=dp), DIMENSION(:), INTENT(IN) :: exc_grad
891 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
892 TYPE(mp_para_env_type), POINTER :: para_env
893
894 CHARACTER(len=1), DIMENSION(3), PARAMETER :: label = ["x", "y", "z"]
895
896 INTEGER :: i, iatom, iw, j
897 REAL(kind=dp), DIMENSION(3) :: center, displacement, grad, grad_sum
898 REAL(kind=dp), DIMENSION(3, 3) :: molecular_virial
899
900 cpassert(ASSOCIATED(particle_set))
901 cpassert(SIZE(exc_grad) == 3*SIZE(particle_set))
902
903 IF (para_env%mepos /= 0) RETURN
904
905 center = 0.0_dp
906 DO iatom = 1, SIZE(particle_set)
907 center = center + particle_set(iatom)%r
908 END DO
909 center = center/real(SIZE(particle_set), dp)
910
911 grad_sum = 0.0_dp
912 molecular_virial = 0.0_dp
913 DO iatom = 1, SIZE(particle_set)
914 grad = exc_grad(3*iatom - 2:3*iatom)
915 displacement = particle_set(iatom)%r - center
916 grad_sum = grad_sum + grad
917 DO i = 1, 3
918 DO j = 1, 3
919 molecular_virial(i, j) = molecular_virial(i, j) + grad(i)*displacement(j)
920 END DO
921 END DO
922 END DO
923
924 iw = cp_logger_get_default_io_unit()
925 IF (iw <= 0) RETURN
926
927 WRITE (unit=iw, fmt="(/,T2,A)") &
928 "GAUXC| Molecular XC gradient virial diagnostic [a.u.]"
929 WRITE (unit=iw, fmt="(T2,A,T20,A,T40,A,T60,A)") "GAUXC|", "x", "y", "z"
930 DO i = 1, 3
931 WRITE (unit=iw, fmt="(T2,A,1X,A1,3(1X,ES19.11))") &
932 "GAUXC|", label(i), molecular_virial(i, :)
933 END DO
934 WRITE (unit=iw, fmt="(T2,A,1X,ES19.11)") &
935 "GAUXC| Molecular XC gradient virial 1/3 Trace", &
936 (molecular_virial(1, 1) + molecular_virial(2, 2) + molecular_virial(3, 3))/3.0_dp
937 WRITE (unit=iw, fmt="(T2,A,3(1X,ES19.11))") &
938 "GAUXC| Molecular XC gradient sum", grad_sum
939 WRITE (unit=iw, fmt="(T2,A)") &
940 "GAUXC| Diagnostic only; this is not an analytical periodic stress tensor."
941
942 END SUBROUTINE print_gauxc_molecular_virial
943
944! **************************************************************************************************
945!> \brief Return information about the Skala functional
946!> \param functional section containing the SKALA subsection
947!> \param lsd if you are using lsd or lda
948!> \param reference the reference to the article where the functional is explained
949!> \param shortform the short definition of the functional
950!> \param needs the flags corresponding to the inputs needed by this
951!> functional are set to true (the flags not needed aren't touched)
952!> \param max_deriv the maximal derivative available
953! **************************************************************************************************
954 SUBROUTINE skala_info(functional, lsd, reference, shortform, needs, max_deriv)
955 TYPE(section_vals_type), POINTER :: functional
956 LOGICAL, INTENT(in) :: lsd
957 CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
958 TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
959 INTEGER, INTENT(out), OPTIONAL :: max_deriv
960
961 CHARACTER(len=default_path_length) :: model_key, model_name
962 CHARACTER(len=default_string_length) :: xc_fun_key, xc_fun_name
963 INTEGER :: gapw_representation
964 LOGICAL :: native_grid
965
966 CALL section_vals_val_get(functional, "FUNCTIONAL", c_val=xc_fun_name)
967 CALL section_vals_val_get(functional, "MODEL", c_val=model_name)
968 CALL section_vals_val_get(functional, "NATIVE_GRID", l_val=native_grid)
969 CALL section_vals_val_get(functional, "PSEUDOPOTENTIAL_GAPW_REPRESENTATION", &
970 i_val=gapw_representation)
971 native_grid = native_grid .OR. gapw_representation == skala_gapw_paw_one_center
972 model_key = adjustl(model_name)
973 xc_fun_key = adjustl(xc_fun_name)
974 CALL uppercase(model_key)
975 CALL uppercase(xc_fun_key)
976
977 IF (PRESENT(reference)) THEN
978 IF (trim(model_key) == "NONE" .OR. trim(model_key) == "" .OR. &
979 trim(model_key) == trim(xc_fun_key)) THEN
980 reference = "Functional computed by GauXC (underlying: "//trim(xc_fun_name)//")"
981 ELSE
982 reference = "Functional computed by GauXC Skala model "//trim(model_name)
983 END IF
984 END IF
985 IF (PRESENT(shortform)) THEN
986 IF (trim(model_key) == "NONE" .OR. trim(model_key) == "" .OR. &
987 trim(model_key) == trim(xc_fun_key)) THEN
988 shortform = "GAUXC ("//trim(xc_fun_name)//")"
989 ELSE
990 shortform = "GAUXC Skala"
991 END IF
992 END IF
993 IF (PRESENT(needs)) THEN
994 IF (native_grid .AND. trim(model_key) /= "NONE" .AND. trim(model_key) /= "" .AND. &
995 trim(model_key) /= trim(xc_fun_key)) THEN
996 IF (lsd) THEN
997 needs%rho_spin = .true.
998 needs%drho_spin = .true.
999 needs%tau_spin = .true.
1000 ELSE
1001 needs%rho = .true.
1002 needs%drho = .true.
1003 needs%tau = .true.
1004 END IF
1005 ELSE
1006 needs%rho = .true.
1007 IF (lsd) THEN
1008 needs%rho_spin = .true.
1009 END IF
1010 END IF
1011 END IF
1012 IF (PRESENT(max_deriv)) max_deriv = 1
1013
1014 END SUBROUTINE skala_info
1015
1016! GauXC uses replicated dense density and VXC matrices. The DBCSR density matrix
1017! is distributed over MPI ranks, so apply_gauxc allreduces the dense copy before
1018! passing it to GauXC.
1019
1020! **************************************************************************************************
1021!> \brief ...
1022!> \param qs_env ...
1023!> \param xc_section ...
1024!> \param calculate_forces ...
1025! **************************************************************************************************
1026 SUBROUTINE apply_gauxc(qs_env, xc_section, calculate_forces)
1027 TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
1028 TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
1029 LOGICAL, INTENT(IN) :: calculate_forces
1030
1031 CHARACTER(len=*), PARAMETER :: nonlocal_vdw_abort_message = &
1032 "GauXC does not support non-local VDW_POTENTIAL corrections. "// &
1033 "Use an additive PAIR_POTENTIAL dispersion correction or disable GauXC."
1034 REAL(kind=dp), PARAMETER :: gapw_fd_gradient_dx = 1.0e-4_dp
1035
1036 CHARACTER(len=default_path_length) :: model_key, model_name, output_path
1037 CHARACTER(len=default_string_length) :: gradient_runtime, gradient_runtime_key, grid_key, &
1038 int_exec_space_key, lwd_kernel_key, pruning_key, skala_runtime, skala_runtime_key, &
1039 xc_fun_key
1040 INTEGER :: atom_chunk_size, env_status, img, ispin, &
1041 nimages
1042 LOGICAL :: atom_chunk_size_explicit, do_kpoints, gapw_method, gapw_paw_pseudopotentials, &
1043 gapw_pseudopotentials, grid_explicit, hdf5_output, is_periodic, molecular_virial, &
1044 molecular_virial_debug, need_xc_gradient, periodic_reference, pruning_explicit, &
1045 use_skala_model, write_hdf5_output
1046 REAL(kind=dp) :: molecular_virial_debug_dx
1047 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: density_scalar, density_zeta
1048 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1049 TYPE(cell_type), POINTER :: cell
1050 TYPE(cp_gauxc_cache_params) :: params
1051 TYPE(cp_gauxc_cache_type), POINTER :: cache
1052 TYPE(cp_gauxc_status_type) :: gauxc_status
1053 TYPE(cp_gauxc_xc_gradient_type) :: exc_grad
1054 TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
1055 TYPE(dbcsr_p_type) :: vxc_zeta_tmp
1056 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
1057 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
1058 TYPE(dft_control_type), POINTER :: dft_control
1059 TYPE(mp_para_env_type), POINTER :: para_env
1060 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1061 TYPE(qs_energy_type), POINTER :: energy
1062 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1063 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1064 TYPE(qs_ks_env_type), POINTER :: ks_env
1065 TYPE(qs_rho_type), POINTER :: rho, rho_use, rho_xc
1066 TYPE(qs_scf_env_type), POINTER :: scf_env
1067 TYPE(section_vals_type), POINTER :: gauxc_functional_section
1068
1069 NULLIFY ( &
1070 atomic_kind_set, &
1071 cell, &
1072 dft_control, &
1073 energy, &
1074 force, &
1075 ks_env, &
1076 matrix_vxc, &
1077 para_env, &
1078 particle_set, &
1079 qs_kind_set, &
1080 rho, &
1081 rho_use, &
1082 rho_xc, &
1083 rho_ao, &
1084 scf_env)
1085
1086 CALL get_qs_env( &
1087 qs_env, &
1088 cell=cell, &
1089 dft_control=dft_control, &
1090 do_kpoints=do_kpoints, &
1091 energy=energy, &
1092 ks_env=ks_env, &
1093 matrix_vxc=matrix_vxc, &
1094 natom=params%natom, &
1095 atomic_kind_set=atomic_kind_set, &
1096 force=force, &
1097 para_env=para_env, &
1098 particle_set=particle_set, &
1099 qs_kind_set=qs_kind_set, &
1100 rho=rho, &
1101 rho_xc=rho_xc, &
1102 scf_env=scf_env)
1103
1104 gapw_method = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1105 gapw_pseudopotentials = gapw_method .AND. &
1106 gauxc_gapw_has_pseudopotentials(qs_kind_set)
1107 gapw_paw_pseudopotentials = gapw_method .AND. &
1108 gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
1109 IF (dft_control%qs_control%gapw_xc) THEN
1110 cpassert(ASSOCIATED(rho_xc))
1111 rho_use => rho_xc
1112 ELSE
1113 cpassert(ASSOCIATED(rho))
1114 rho_use => rho
1115 END IF
1116 CALL qs_rho_get( &
1117 rho_use, &
1118 rho_ao_kp=rho_ao)
1119
1120 nimages = dft_control%nimages
1121 params%nspins = dft_control%nspins
1122 is_periodic = .false.
1123 IF (ASSOCIATED(cell)) is_periodic = any(cell%perd /= 0)
1124
1125 IF (ASSOCIATED(qs_env%dispersion_env)) THEN
1126 IF (qs_env%dispersion_env%type == xc_vdw_fun_nonloc) THEN
1127 cpabort(nonlocal_vdw_abort_message)
1128 END IF
1129 END IF
1130 NULLIFY (vxc_zeta_tmp%matrix)
1131
1132 gauxc_functional_section => get_gauxc_functional(xc_section)
1133 CALL section_vals_val_get( &
1134 gauxc_functional_section, &
1135 "FUNCTIONAL", &
1136 c_val=params%xc_fun_name)
1137 CALL section_vals_val_get( &
1138 gauxc_functional_section, &
1139 "MODEL", &
1140 c_val=model_name)
1141 CALL section_vals_val_get( &
1142 gauxc_functional_section, &
1143 "GRID", &
1144 c_val=params%grid_type, &
1145 explicit=grid_explicit)
1146 CALL section_vals_val_get( &
1147 gauxc_functional_section, &
1148 "RADIAL_QUADRATURE", &
1149 c_val=params%radial_quadrature)
1150 CALL section_vals_val_get( &
1151 gauxc_functional_section, &
1152 "PRUNING_SCHEME", &
1153 c_val=params%pruning_scheme, &
1154 explicit=pruning_explicit)
1155 CALL section_vals_val_get( &
1156 gauxc_functional_section, &
1157 "BATCH_SIZE", &
1158 i_val=params%batch_size)
1159 CALL section_vals_val_get( &
1160 gauxc_functional_section, &
1161 "DEVICE_RUNTIME_FILL_FRACTION", &
1162 r_val=params%device_runtime_fill_fraction)
1163 CALL section_vals_val_get( &
1164 gauxc_functional_section, &
1165 "MODEL_ATOM_CHUNK_SIZE", &
1166 i_val=atom_chunk_size, &
1167 explicit=atom_chunk_size_explicit)
1168 CALL section_vals_val_get( &
1169 gauxc_functional_section, &
1170 "PERIODIC_REFERENCE", &
1171 l_val=periodic_reference)
1172 CALL section_vals_val_get( &
1173 gauxc_functional_section, &
1174 "MOLECULAR_VIRIAL", &
1175 l_val=molecular_virial)
1176 CALL section_vals_val_get( &
1177 gauxc_functional_section, &
1178 "MOLECULAR_VIRIAL_DEBUG", &
1179 l_val=molecular_virial_debug)
1180 CALL section_vals_val_get( &
1181 gauxc_functional_section, &
1182 "MOLECULAR_VIRIAL_DEBUG_DX", &
1183 r_val=molecular_virial_debug_dx)
1184 CALL section_vals_val_get( &
1185 gauxc_functional_section, &
1186 "LB_EXECUTION_SPACE", &
1187 c_val=params%lb_exec_space)
1188 CALL section_vals_val_get( &
1189 gauxc_functional_section, &
1190 "INT_EXECUTION_SPACE", &
1191 c_val=params%int_exec_space)
1192 CALL section_vals_val_get( &
1193 gauxc_functional_section, &
1194 "LWD_KERNEL", &
1195 c_val=params%lwd_kernel)
1196 CALL section_vals_val_get( &
1197 gauxc_functional_section, &
1198 "SKALA_RUNTIME", &
1199 c_val=skala_runtime)
1200 CALL section_vals_val_get( &
1201 gauxc_functional_section, &
1202 "MODEL_GRADIENT_RUNTIME", &
1203 c_val=gradient_runtime)
1204 CALL section_vals_val_get( &
1205 gauxc_functional_section, &
1206 "OUTPUT_PATH", &
1207 c_val=output_path)
1208
1209 model_key = adjustl(model_name)
1210 CALL uppercase(model_key)
1211 xc_fun_key = adjustl(params%xc_fun_name)
1212 CALL uppercase(xc_fun_key)
1213 skala_runtime_key = adjustl(skala_runtime)
1214 CALL uppercase(skala_runtime_key)
1215 gradient_runtime_key = adjustl(gradient_runtime)
1216 CALL uppercase(gradient_runtime_key)
1217 int_exec_space_key = adjustl(params%int_exec_space)
1218 CALL uppercase(int_exec_space_key)
1219 params%use_gauxc_model = (trim(model_key) /= "" .AND. trim(model_key) /= "NONE" .AND. &
1220 trim(model_key) /= trim(xc_fun_key))
1221 use_skala_model = (index(trim(model_key), "SKALA") > 0)
1222 lwd_kernel_key = adjustl(params%lwd_kernel)
1223 CALL uppercase(lwd_kernel_key)
1224 IF (trim(lwd_kernel_key) == "AUTO") THEN
1225 params%lwd_kernel = "DEFAULT"
1226#if defined(GAUXC_HAS_CUTLASS)
1227 IF (use_skala_model .AND. trim(params%int_exec_space) == "DEVICE") THEN
1228 params%lwd_kernel = "SCHEME1-CUTLASS"
1229 END IF
1230#endif
1231 END IF
1232 params%model_eval_name = model_name
1233 IF (.NOT. params%use_gauxc_model) THEN
1234 ! MODEL NONE and MODEL equal to FUNCTIONAL select conventional GauXC.
1235 params%model_eval_name = "NONE"
1236 END IF
1237 IF (gapw_pseudopotentials .AND. params%use_gauxc_model .AND. .NOT. dft_control%qs_control%gapw_xc .AND. &
1238 .NOT. gapw_paw_pseudopotentials .AND. para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
1239 IF (scf_env%iter_count == 1) THEN
1240 CALL cp_warn( &
1241 __location__, &
1242 "GauXC Skala with METHOD GAPW and GPW_TYPE pseudopotentials evaluates "// &
1243 "the XC term directly on the molecular AO/valence density; no GAPW one-center "// &
1244 "XC correction is used for those regular-grid kinds.")
1245 END IF
1246 END IF
1247 IF (params%device_runtime_fill_fraction <= 0.0_dp .OR. params%device_runtime_fill_fraction > 1.0_dp) THEN
1248 CALL cp_abort(__location__, &
1249 "GAUXC%DEVICE_RUNTIME_FILL_FRACTION must be > 0 and <= 1.")
1250 END IF
1251 IF (atom_chunk_size < -1) THEN
1252 CALL cp_abort(__location__, &
1253 "GAUXC%MODEL_ATOM_CHUNK_SIZE must be -1, zero, or positive.")
1254 END IF
1255 IF (molecular_virial_debug) THEN
1256 IF (molecular_virial_debug_dx <= 0.0_dp) THEN
1257 CALL cp_abort(__location__, &
1258 "GauXC MOLECULAR_VIRIAL_DEBUG_DX must be positive.")
1259 END IF
1260 molecular_virial = .true.
1261 END IF
1262 need_xc_gradient = calculate_forces .OR. molecular_virial
1263 CALL ensure_gauxc_periodic_reference_scope( &
1264 dft_control, cell, qs_kind_set, do_kpoints, periodic_reference)
1265 IF (is_periodic .AND. periodic_reference .AND. para_env%mepos == 0) THEN
1266 IF (ASSOCIATED(scf_env)) THEN
1267 IF (scf_env%iter_count == 1) THEN
1268 CALL cp_warn( &
1269 __location__, &
1270 "GAUXC%PERIODIC_REFERENCE uses GauXC molecular quadrature for isolated validation "// &
1271 "cells. Compact periodic materials require a dedicated periodic GauXC interface.")
1272 END IF
1273 END IF
1274 END IF
1275 IF (params%use_gauxc_model) THEN
1276 IF (has_nlcc(qs_kind_set)) THEN
1277 CALL cp_abort(__location__, &
1278 "GauXC Skala with NLCC pseudopotentials is not implemented. "// &
1279 "The frozen core density would need a SKALA-consistent feature definition.")
1280 END IF
1281 END IF
1282 IF (params%use_gauxc_model) THEN
1283 CALL set_gauxc_model_atom_chunk_env( &
1284 atom_chunk_size, atom_chunk_size_explicit)
1285 IF (.NOT. grid_explicit) params%grid_type = "SUPERFINE"
1286 IF (.NOT. pruning_explicit) params%pruning_scheme = "UNPRUNED"
1287
1288 grid_key = adjustl(params%grid_type)
1289 pruning_key = adjustl(params%pruning_scheme)
1290 CALL uppercase(grid_key)
1291 CALL uppercase(pruning_key)
1292 IF (use_skala_model .AND. need_xc_gradient .AND. &
1293 (trim(grid_key) /= "SUPERFINE" .OR. trim(pruning_key) /= "UNPRUNED")) THEN
1294 CALL cp_warn( &
1295 __location__, &
1296 "GauXC Skala nuclear gradients are sensitive to the GauXC molecular grid. "// &
1297 "Use GRID SUPERFINE and PRUNING_SCHEME UNPRUNED for quantitative force checks.")
1298 END IF
1299 IF (trim(model_key) == "SKALA") THEN
1300 model_name = ""
1301 env_status = 1
1302 IF (trim(int_exec_space_key) == "DEVICE") THEN
1303 CALL get_environment_variable("GAUXC_SKALA_CUDA_MODEL", model_name, status=env_status)
1304 END IF
1305 IF (env_status /= 0 .OR. len_trim(model_name) == 0) THEN
1306 CALL get_environment_variable("GAUXC_SKALA_MODEL", model_name, status=env_status)
1307 END IF
1308 IF (env_status /= 0 .OR. len_trim(model_name) == 0) THEN
1309 IF (trim(int_exec_space_key) == "DEVICE") THEN
1310 CALL cp_abort( &
1311 __location__, &
1312 "MODEL SKALA with DEVICE execution requires GAUXC_SKALA_CUDA_MODEL or "// &
1313 "GAUXC_SKALA_MODEL")
1314 ELSE
1315 cpabort("MODEL SKALA requires the GAUXC_SKALA_MODEL environment variable")
1316 END IF
1317 END IF
1318 params%model_eval_name = model_name
1319 END IF
1320 END IF
1321 SELECT CASE (trim(skala_runtime_key))
1322 CASE ("AUTO")
1323 params%use_self_runtime = use_skala_model .AND. para_env%num_pe > 1 .AND. params%nspins > 1
1324 CASE ("MPI")
1325 params%use_self_runtime = .false.
1326 CASE ("SELF")
1327 params%use_self_runtime = use_skala_model .AND. para_env%num_pe > 1
1328 CASE DEFAULT
1329 CALL cp_abort(__location__, "Unknown GAUXC%SKALA_RUNTIME value.")
1330 END SELECT
1331 IF (.NOT. use_skala_model) params%use_self_runtime = .false.
1332 SELECT CASE (trim(gradient_runtime_key))
1333 CASE ("AUTO", "SELF")
1334 params%use_gradient_mpi_runtime = .false.
1335 params%use_gradient_self_runtime = need_xc_gradient .AND. params%use_gauxc_model .AND. &
1336 para_env%num_pe > 1 .AND. .NOT. params%use_self_runtime
1337 CASE ("MPI")
1338 params%use_gradient_mpi_runtime = need_xc_gradient .AND. params%use_gauxc_model .AND. para_env%num_pe > 1
1339 params%use_gradient_self_runtime = .false.
1340 CASE DEFAULT
1341 CALL cp_abort(__location__, "Unknown GAUXC%MODEL_GRADIENT_RUNTIME value.")
1342 END SELECT
1343 IF (.NOT. params%use_gauxc_model) THEN
1344 params%use_gradient_mpi_runtime = .false.
1345 params%use_gradient_self_runtime = .false.
1346 END IF
1347 IF (use_skala_model .AND. para_env%num_pe > 1 .AND. .NOT. params%use_self_runtime .AND. &
1348 para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
1349 IF (scf_env%iter_count == 1) THEN
1350 CALL cp_warn( &
1351 __location__, &
1352 "GAUXC%SKALA_RUNTIME uses the MPI communicator for energy/VXC. "// &
1353 "SKALA Torch atom chunks can be distributed across MPI ranks; "// &
1354 "set GAUXC_ONEDFT_DISTRIBUTED_TORCH=0 to force rank-0 Torch inference.")
1355 END IF
1356 END IF
1357
1358 ! After creating the basisset, we will have to check max_l>3 as a further condition
1359 params%use_fd_gradient = gapw_method .AND. need_xc_gradient
1360
1361 IF (.NOT. ASSOCIATED(qs_env%gauxc_cache)) ALLOCATE (qs_env%gauxc_cache)
1362 cache => qs_env%gauxc_cache
1363 CALL gauxc_cache_init( &
1364 cache, &
1365 params, &
1366 para_env, &
1367 particle_set, &
1368 qs_kind_set, &
1369 gauxc_status)
1370
1371 hdf5_output = (trim(output_path) /= "")
1372 write_hdf5_output = hdf5_output .AND. para_env%mepos == 0
1373 IF (write_hdf5_output .AND. ASSOCIATED(scf_env)) THEN
1374 write_hdf5_output = scf_env%iter_count == 1
1375 END IF
1376 IF (write_hdf5_output) THEN
1377 CALL gauxc_write_molecule_hdf5( &
1378 cache%molecule, &
1379 output_path, &
1380 "molecule.h5", &
1381 "molecule", &
1382 gauxc_status)
1383 CALL gauxc_check_status(gauxc_status)
1384 CALL gauxc_write_basisset_hdf5( &
1385 cache%basisset, &
1386 output_path, &
1387 "basisset.h5", &
1388 "basisset", &
1389 gauxc_status)
1390 CALL gauxc_check_status(gauxc_status)
1391 END IF
1392
1393 IF (qs_env%run_rtp) THEN
1394 cpabort("GAUXC XC energy currently does not support real-time propagation")
1395 END IF
1396
1397 energy%exc = 0
1398
1399 IF (ASSOCIATED(matrix_vxc)) CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1400 CALL dbcsr_allocate_matrix_set(matrix_vxc, params%nspins)
1401
1402 DO img = 1, nimages
1403 IF (img > 1) THEN
1404 cpabort("UNIMPLEMENTED: Handling nimg>1 in k-point integration")
1405 END IF
1406 CALL dbcsr_to_dense(rho_ao(1, img), density_scalar, para_env)
1407 CALL para_env%sum(density_scalar)
1408 IF (params%nspins == 1) THEN
1409 gauxc_xc_result = gauxc_compute_xc( &
1410 cache%integrator, &
1411 density_scalar, &
1412 nspins=params%nspins, &
1413 status=gauxc_status, &
1414 model=trim(params%model_eval_name))
1415 CALL gauxc_check_status(gauxc_status)
1416 IF (need_xc_gradient) THEN
1417 IF (params%use_fd_gradient) THEN
1418 CALL gauxc_xc_gradient_fd( &
1419 particle_set, qs_kind_set, density_scalar, params%nspins, params%model_eval_name, &
1420 params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1421 params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1422 params%device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
1423 exc_grad%exc_grad)
1424 ELSE IF (params%use_gradient_self_runtime) THEN
1425 exc_grad = gauxc_compute_xc_gradient( &
1426 cache%gradient_integrator, &
1427 density_scalar, &
1428 nspins=params%nspins, &
1429 natom=params%natom, &
1430 status=gauxc_status, &
1431 model=trim(params%model_eval_name))
1432 ELSE
1433 exc_grad = gauxc_compute_xc_gradient( &
1434 cache%integrator, &
1435 density_scalar, &
1436 nspins=params%nspins, &
1437 natom=params%natom, &
1438 status=gauxc_status, &
1439 model=trim(params%model_eval_name))
1440 END IF
1441 CALL gauxc_check_status(gauxc_status)
1442 IF (calculate_forces) THEN
1443 CALL add_gauxc_gradient_to_force( &
1444 exc_grad%exc_grad, &
1445 force, &
1446 atomic_kind_set, &
1447 para_env)
1448 END IF
1449 IF (molecular_virial) THEN
1450 CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1451 END IF
1452 IF (molecular_virial_debug) THEN
1453 CALL debug_gauxc_molecular_virial( &
1454 exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, params%nspins, &
1455 params%model_eval_name, params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1456 params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1457 params%device_runtime_fill_fraction, molecular_virial_debug_dx, para_env)
1458 END IF
1459 DEALLOCATE (exc_grad%exc_grad)
1460 END IF
1461 ELSE
1462 cpassert(params%nspins == 2)
1463 ! In here:
1464 ! scalar <- rho_ao(1, :) + rho_ao(2, :)
1465 ! zeta <- rho_ao(1, :) - rho_ao(2, :)
1466 CALL dbcsr_to_dense(rho_ao(2, img), density_zeta, para_env)
1467 CALL para_env%sum(density_zeta)
1468 ! Do NOT reorder the following lines!
1469 density_scalar(:, :) = density_scalar(:, :) + density_zeta(:, :)
1470 ! Factor two because the next line is evaluated after the above line.
1471 ! We need to subtract density_zeta once to undo the above line and
1472 ! a second time because that is what UKS requires.
1473 ! This style lowers memory footprint.
1474 density_zeta(:, :) = density_scalar(:, :) - 2.0_dp*density_zeta(:, :)
1475 gauxc_xc_result = gauxc_compute_xc( &
1476 cache%integrator, &
1477 density_scalar, &
1478 density_zeta, &
1479 params%nspins, &
1480 gauxc_status, &
1481 model=trim(params%model_eval_name))
1482 CALL gauxc_check_status(gauxc_status)
1483 IF (need_xc_gradient) THEN
1484 IF (params%use_fd_gradient) THEN
1485 CALL gauxc_xc_gradient_fd( &
1486 particle_set, qs_kind_set, density_scalar, params%nspins, params%model_eval_name, &
1487 params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1488 params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1489 params%device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
1490 exc_grad%exc_grad, density_zeta=density_zeta)
1491 ELSE IF (params%use_gradient_self_runtime) THEN
1492 exc_grad = gauxc_compute_xc_gradient( &
1493 cache%gradient_integrator, &
1494 density_scalar, &
1495 density_zeta, &
1496 params%nspins, &
1497 params%natom, &
1498 gauxc_status, &
1499 model=trim(params%model_eval_name))
1500 ELSE
1501 exc_grad = gauxc_compute_xc_gradient( &
1502 cache%integrator, &
1503 density_scalar, &
1504 density_zeta, &
1505 params%nspins, &
1506 params%natom, &
1507 gauxc_status, &
1508 model=trim(params%model_eval_name))
1509 END IF
1510 CALL gauxc_check_status(gauxc_status)
1511 IF (calculate_forces) THEN
1512 CALL add_gauxc_gradient_to_force( &
1513 exc_grad%exc_grad, &
1514 force, &
1515 atomic_kind_set, &
1516 para_env)
1517 END IF
1518 IF (molecular_virial) THEN
1519 CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1520 END IF
1521 IF (molecular_virial_debug) THEN
1522 CALL debug_gauxc_molecular_virial( &
1523 exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, params%nspins, &
1524 params%model_eval_name, params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1525 params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1526 params%device_runtime_fill_fraction, molecular_virial_debug_dx, para_env, &
1527 density_zeta=density_zeta)
1528 END IF
1529 DEALLOCATE (exc_grad%exc_grad)
1530 END IF
1531 END IF
1532
1533 energy%exc = energy%exc + gauxc_xc_result%exc
1534
1535 IF (params%nspins == 1) THEN
1536 IF (img == 1) THEN
1537 matrix_vxc(1) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(1, img))
1538 ELSE
1539 cpabort("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1540 END IF
1541 ELSE
1542 cpassert(params%nspins == 2)
1543 ! Transform derivatives from total/spin density back to alpha/beta channels.
1544 vxc_zeta_tmp = dense_to_dbcsr(gauxc_xc_result%vxc_zeta, rho_ao(1, img))
1545 IF (img == 1) THEN
1546 DO ispin = 1, 2
1547 matrix_vxc(ispin) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(ispin, 1))
1548 CALL dbcsr_add( &
1549 matrix_vxc(ispin)%matrix, &
1550 vxc_zeta_tmp%matrix, &
1551 1.0_dp, &
1552 ! 1.0 for ispin==1, -1.0 for ispin==2
1553 1.0_dp - real(ispin - 1, dp)*2.0_dp)
1554 END DO
1555 ELSE
1556 cpabort("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1557 END IF
1558 CALL dbcsr_release(vxc_zeta_tmp%matrix)
1559 DEALLOCATE (vxc_zeta_tmp%matrix)
1560 END IF
1561 END DO
1562
1563 DEALLOCATE (density_scalar)
1564 IF (ALLOCATED(density_zeta)) DEALLOCATE (density_zeta)
1565 DEALLOCATE (gauxc_xc_result%vxc_scalar)
1566 IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
1567
1568 CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
1569 DO ispin = 1, params%nspins
1570 CALL dbcsr_finalize(matrix_vxc(ispin)%matrix)
1571 END DO
1572
1573 END SUBROUTINE apply_gauxc
1574
1575END MODULE xc_gauxc_functional
cp_dbcsr_api::dbcsr_create
Definition cp_dbcsr_api.F:194
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:77
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:85
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
cp_dbcsr_api::dbcsr_get_readonly_block_p
subroutine, public dbcsr_get_readonly_block_p(matrix, row, col, block, found, row_size, col_size)
Like dbcsr_get_block_p() but with matrix being INTENT(IN). When invoking this routine,...
Definition cp_dbcsr_api.F:718
cp_dbcsr_api::dbcsr_get_matrix_type
character function, public dbcsr_get_matrix_type(matrix)
...
Definition cp_dbcsr_api.F:861
cp_dbcsr_api::dbcsr_get_info
subroutine, public dbcsr_get_info(matrix, nblkrows_total, nblkcols_total, nfullrows_total, nfullcols_total, nblkrows_local, nblkcols_local, nfullrows_local, nfullcols_local, my_prow, my_pcol, local_rows, local_cols, proc_row_dist, proc_col_dist, row_blk_size, col_blk_size, row_blk_offset, col_blk_offset, distribution, name, matrix_type, group)
...
Definition cp_dbcsr_api.F:808
cp_dbcsr_api::dbcsr_get_stored_coordinates
subroutine, public dbcsr_get_stored_coordinates(matrix, row, column, processor)
...
Definition cp_dbcsr_api.F:911
cp_dbcsr_api::dbcsr_init_p
subroutine, public dbcsr_init_p(matrix)
...
Definition cp_dbcsr_api.F:207
cp_dbcsr_api::dbcsr_work_create
subroutine, public dbcsr_work_create(matrix, nblks_guess, sizedata_guess, n, work_mutable)
...
Definition cp_dbcsr_api.F:1307
cp_dbcsr_api::dbcsr_finalize
subroutine, public dbcsr_finalize(matrix)
...
Definition cp_dbcsr_api.F:673
cp_dbcsr_api::dbcsr_release
subroutine, public dbcsr_release(matrix)
...
Definition cp_dbcsr_api.F:1142
cp_dbcsr_api::dbcsr_put_block
subroutine, public dbcsr_put_block(matrix, row, col, block, summation)
...
Definition cp_dbcsr_api.F:1125
cp_dbcsr_api::dbcsr_add
subroutine, public dbcsr_add(matrix_a, matrix_b, alpha_scalar, beta_scalar)
...
Definition cp_dbcsr_api.F:254
cp_dbcsr_api::dbcsr_distribution_get
subroutine, public dbcsr_distribution_get(dist, row_dist, col_dist, nrows, ncols, has_threads, group, mynode, numnodes, nprows, npcols, myprow, mypcol, pgrid, subgroups_defined, prow_group, pcol_group)
...
Definition cp_dbcsr_api.F:572
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::xc_vdw_fun_nonloc
integer, parameter, public xc_vdw_fun_nonloc
Definition input_constants.F:621
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals2
type(section_vals_type) function, pointer, public section_vals_get_subs_vals2(section_vals, i_section, i_rep_section)
returns the values of the n-th non default subsection (null if no such section exists (not so many no...
Definition input_section_types.F:791
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
message_passing::mp_comm_self
type(mp_comm_type), parameter, public mp_comm_self
Definition message_passing.F:816
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_energy_types
Definition qs_energy_types.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:539
qs_force_types
Definition qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::has_nlcc
logical function, public has_nlcc(qs_kind_set)
finds if a given qs run needs to use nlcc
Definition qs_kind_types.F:4286
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, exc_accint, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:591
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
xc_gauxc_cache
Definition xc_gauxc_cache.F:12
xc_gauxc_cache::gauxc_cache_init
subroutine, public gauxc_cache_init(cache, params, para_env, particle_set, qs_kind_set, status)
...
Definition xc_gauxc_cache.F:116
xc_gauxc_functional
Definition xc_gauxc_functional.F:12
xc_gauxc_functional::gauxc_gapw_has_paw_pseudopotentials
logical function, public gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
Return whether GauXC GAPW mode sees pseudopotential one-center GAPW kinds.
Definition xc_gauxc_functional.F:450
xc_gauxc_functional::xc_section_uses_gauxc
logical function, public xc_section_uses_gauxc(xc_section)
...
Definition xc_gauxc_functional.F:361
xc_gauxc_functional::gauxc_gapw_has_pseudopotentials
logical function, public gauxc_gapw_has_pseudopotentials(qs_kind_set)
Return whether GauXC GAPW mode sees pseudopotential kinds.
Definition xc_gauxc_functional.F:392
xc_gauxc_functional::skala_info
subroutine, public skala_info(functional, lsd, reference, shortform, needs, max_deriv)
Return information about the Skala functional.
Definition xc_gauxc_functional.F:955
xc_gauxc_functional::gauxc_gapw_all_all_electron
logical function, public gauxc_gapw_all_all_electron(qs_kind_set)
Return whether every GAPW kind uses an all-electron potential.
Definition xc_gauxc_functional.F:103
xc_gauxc_functional::gauxc_gapw_all_pseudopotentials
logical function, public gauxc_gapw_all_pseudopotentials(qs_kind_set)
Return whether every GAPW kind uses a pseudopotential.
Definition xc_gauxc_functional.F:421
xc_gauxc_functional::apply_gauxc
subroutine, public apply_gauxc(qs_env, xc_section, calculate_forces)
...
Definition xc_gauxc_functional.F:1027
xc_input_constants
input constants for xc
Definition xc_input_constants.F:10
xc_input_constants::skala_gapw_paw_one_center
integer, parameter, public skala_gapw_paw_one_center
Definition xc_input_constants.F:29
xc_rho_cflags_types
contains the structure
Definition xc_rho_cflags_types.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:758
cp_dbcsr_api::dbcsr_distribution_type
Definition cp_dbcsr_api.F:182
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
external_potential_types::all_potential_type
Definition external_potential_types.F:70
external_potential_types::gth_potential_type
Definition external_potential_types.F:116
external_potential_types::sgp_potential_type
Definition external_potential_types.F:174
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:767
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:137
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:99
xc_gauxc_cache::cp_gauxc_cache_params
Definition xc_gauxc_cache.F:31
xc_gauxc_cache::cp_gauxc_cache_type
Definition xc_gauxc_cache.F:52
xc_rho_cflags_types::xc_rho_cflags_type
contains a flag for each component of xc_rho_set, so that you can use it to tell which components you...
Definition xc_rho_cflags_types.F:48