(git:08825b9)
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 Evaluate molecular GauXC, optionally for an explicit configuration density.
1022!> \param qs_env ...
1023!> \param xc_section ...
1024!> \param calculate_forces ...
1025!> \param rho_ao_external Optional spin/image AO density; supply all three external arguments together
1026!> \param exc_external Configuration XC energy, leaving the baseline energy unchanged
1027!> \param matrix_vxc_external Initially null; caller owns the returned spin AO potential matrices
1028!> \note External-density mode preserves baseline energy and KS storage and shares
1029!> the geometry/model cache.
1030! **************************************************************************************************
1031 SUBROUTINE apply_gauxc(qs_env, xc_section, calculate_forces, &
1032 rho_ao_external, exc_external, matrix_vxc_external)
1033 TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
1034 TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
1035 LOGICAL, INTENT(IN) :: calculate_forces
1036 TYPE(dbcsr_p_type), DIMENSION(:, :), INTENT(IN), &
1037 OPTIONAL, POINTER :: rho_ao_external
1038 REAL(kind=dp), INTENT(OUT), OPTIONAL :: exc_external
1039 TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT), &
1040 OPTIONAL, POINTER :: matrix_vxc_external
1041
1042 CHARACTER(len=*), PARAMETER :: nonlocal_vdw_abort_message = &
1043 "GauXC does not support non-local VDW_POTENTIAL corrections. "// &
1044 "Use an additive PAIR_POTENTIAL dispersion correction or disable GauXC."
1045 REAL(kind=dp), PARAMETER :: gapw_fd_gradient_dx = 1.0e-4_dp
1046
1047 CHARACTER(len=default_path_length) :: model_key, model_name, output_path
1048 CHARACTER(len=default_string_length) :: gradient_runtime, gradient_runtime_key, grid_key, &
1049 int_exec_space_key, lwd_kernel_key, pruning_key, skala_runtime, skala_runtime_key, &
1050 xc_fun_key
1051 INTEGER :: atom_chunk_size, env_status, img, ispin, &
1052 nimages
1053 LOGICAL :: atom_chunk_size_explicit, do_kpoints, external_density, gapw_method, &
1054 gapw_paw_pseudopotentials, gapw_pseudopotentials, grid_explicit, hdf5_output, &
1055 is_periodic, low_spin_scope, molecular_virial, molecular_virial_debug, need_xc_gradient, &
1056 periodic_reference, pruning_explicit, use_skala_model, write_hdf5_output
1057 REAL(kind=dp) :: exc_value, molecular_virial_debug_dx
1058 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: density_scalar, density_zeta
1059 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1060 TYPE(cell_type), POINTER :: cell
1061 TYPE(cp_gauxc_cache_params) :: params
1062 TYPE(cp_gauxc_cache_type), POINTER :: cache
1063 TYPE(cp_gauxc_status_type) :: gauxc_status
1064 TYPE(cp_gauxc_xc_gradient_type) :: exc_grad
1065 TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
1066 TYPE(dbcsr_p_type) :: vxc_zeta_tmp
1067 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
1068 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
1069 TYPE(dft_control_type), POINTER :: dft_control
1070 TYPE(mp_para_env_type), POINTER :: para_env
1071 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1072 TYPE(qs_energy_type), POINTER :: energy
1073 TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1074 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1075 TYPE(qs_ks_env_type), POINTER :: ks_env
1076 TYPE(qs_rho_type), POINTER :: rho, rho_use, rho_xc
1077 TYPE(qs_scf_env_type), POINTER :: scf_env
1078 TYPE(section_vals_type), POINTER :: gauxc_functional_section
1079
1080 external_density = PRESENT(rho_ao_external)
1081 cpassert(PRESENT(exc_external) .EQV. external_density)
1082 cpassert(PRESENT(matrix_vxc_external) .EQV. external_density)
1083 IF (external_density) THEN
1084 cpassert(ASSOCIATED(rho_ao_external))
1085 cpassert(.NOT. ASSOCIATED(matrix_vxc_external))
1086 END IF
1087
1088 NULLIFY ( &
1089 atomic_kind_set, &
1090 cell, &
1091 dft_control, &
1092 energy, &
1093 force, &
1094 ks_env, &
1095 matrix_vxc, &
1096 para_env, &
1097 particle_set, &
1098 qs_kind_set, &
1099 rho, &
1100 rho_use, &
1101 rho_xc, &
1102 rho_ao, &
1103 scf_env)
1104
1105 CALL get_qs_env( &
1106 qs_env, &
1107 cell=cell, &
1108 dft_control=dft_control, &
1109 do_kpoints=do_kpoints, &
1110 energy=energy, &
1111 ks_env=ks_env, &
1112 matrix_vxc=matrix_vxc, &
1113 natom=params%natom, &
1114 atomic_kind_set=atomic_kind_set, &
1115 force=force, &
1116 para_env=para_env, &
1117 particle_set=particle_set, &
1118 qs_kind_set=qs_kind_set, &
1119 rho=rho, &
1120 rho_xc=rho_xc, &
1121 scf_env=scf_env)
1122
1123 gapw_method = dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc
1124 gapw_pseudopotentials = gapw_method .AND. &
1125 gauxc_gapw_has_pseudopotentials(qs_kind_set)
1126 gapw_paw_pseudopotentials = gapw_method .AND. &
1127 gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
1128 IF (dft_control%qs_control%gapw_xc) THEN
1129 cpassert(ASSOCIATED(rho_xc))
1130 rho_use => rho_xc
1131 ELSE
1132 cpassert(ASSOCIATED(rho))
1133 rho_use => rho
1134 END IF
1135 CALL qs_rho_get( &
1136 rho_use, &
1137 rho_ao_kp=rho_ao)
1138
1139 IF (external_density) THEN
1140 rho_ao => rho_ao_external
1141 ! get_qs_env returned the baseline pointer. Never release or replace it.
1142 NULLIFY (matrix_vxc)
1143 END IF
1144 nimages = dft_control%nimages
1145 params%nspins = dft_control%nspins
1146 IF (external_density) THEN
1147 cpassert(SIZE(rho_ao, 1) == params%nspins)
1148 cpassert(SIZE(rho_ao, 2) == nimages)
1149 END IF
1150 low_spin_scope = external_density .OR. (dft_control%low_spin_roks .AND. gapw_method)
1151 is_periodic = .false.
1152 IF (ASSOCIATED(cell)) is_periodic = any(cell%perd /= 0)
1153
1154 IF (ASSOCIATED(qs_env%dispersion_env)) THEN
1155 IF (qs_env%dispersion_env%type == xc_vdw_fun_nonloc) THEN
1156 cpabort(nonlocal_vdw_abort_message)
1157 END IF
1158 END IF
1159 NULLIFY (vxc_zeta_tmp%matrix)
1160
1161 gauxc_functional_section => get_gauxc_functional(xc_section)
1162 CALL section_vals_val_get( &
1163 gauxc_functional_section, &
1164 "FUNCTIONAL", &
1165 c_val=params%xc_fun_name)
1166 CALL section_vals_val_get( &
1167 gauxc_functional_section, &
1168 "MODEL", &
1169 c_val=model_name)
1170 CALL section_vals_val_get( &
1171 gauxc_functional_section, &
1172 "GRID", &
1173 c_val=params%grid_type, &
1174 explicit=grid_explicit)
1175 CALL section_vals_val_get( &
1176 gauxc_functional_section, &
1177 "RADIAL_QUADRATURE", &
1178 c_val=params%radial_quadrature)
1179 CALL section_vals_val_get( &
1180 gauxc_functional_section, &
1181 "PRUNING_SCHEME", &
1182 c_val=params%pruning_scheme, &
1183 explicit=pruning_explicit)
1184 CALL section_vals_val_get( &
1185 gauxc_functional_section, &
1186 "BATCH_SIZE", &
1187 i_val=params%batch_size)
1188 CALL section_vals_val_get( &
1189 gauxc_functional_section, &
1190 "DEVICE_RUNTIME_FILL_FRACTION", &
1191 r_val=params%device_runtime_fill_fraction)
1192 CALL section_vals_val_get( &
1193 gauxc_functional_section, &
1194 "MODEL_ATOM_CHUNK_SIZE", &
1195 i_val=atom_chunk_size, &
1196 explicit=atom_chunk_size_explicit)
1197 CALL section_vals_val_get( &
1198 gauxc_functional_section, &
1199 "PERIODIC_REFERENCE", &
1200 l_val=periodic_reference)
1201 CALL section_vals_val_get( &
1202 gauxc_functional_section, &
1203 "MOLECULAR_VIRIAL", &
1204 l_val=molecular_virial)
1205 CALL section_vals_val_get( &
1206 gauxc_functional_section, &
1207 "MOLECULAR_VIRIAL_DEBUG", &
1208 l_val=molecular_virial_debug)
1209 CALL section_vals_val_get( &
1210 gauxc_functional_section, &
1211 "MOLECULAR_VIRIAL_DEBUG_DX", &
1212 r_val=molecular_virial_debug_dx)
1213 CALL section_vals_val_get( &
1214 gauxc_functional_section, &
1215 "LB_EXECUTION_SPACE", &
1216 c_val=params%lb_exec_space)
1217 CALL section_vals_val_get( &
1218 gauxc_functional_section, &
1219 "INT_EXECUTION_SPACE", &
1220 c_val=params%int_exec_space)
1221 CALL section_vals_val_get( &
1222 gauxc_functional_section, &
1223 "LWD_KERNEL", &
1224 c_val=params%lwd_kernel)
1225 CALL section_vals_val_get( &
1226 gauxc_functional_section, &
1227 "SKALA_RUNTIME", &
1228 c_val=skala_runtime)
1229 CALL section_vals_val_get( &
1230 gauxc_functional_section, &
1231 "MODEL_GRADIENT_RUNTIME", &
1232 c_val=gradient_runtime)
1233 CALL section_vals_val_get( &
1234 gauxc_functional_section, &
1235 "FORCE_NUMERICAL_GRADIENT", &
1236 l_val=params%force_fd_gradient)
1237 CALL section_vals_val_get( &
1238 gauxc_functional_section, &
1239 "OUTPUT_PATH", &
1240 c_val=output_path)
1241
1242 model_key = adjustl(model_name)
1243 CALL uppercase(model_key)
1244 xc_fun_key = adjustl(params%xc_fun_name)
1245 CALL uppercase(xc_fun_key)
1246 skala_runtime_key = adjustl(skala_runtime)
1247 CALL uppercase(skala_runtime_key)
1248 gradient_runtime_key = adjustl(gradient_runtime)
1249 CALL uppercase(gradient_runtime_key)
1250 int_exec_space_key = adjustl(params%int_exec_space)
1251 CALL uppercase(int_exec_space_key)
1252 params%use_gauxc_model = (trim(model_key) /= "" .AND. trim(model_key) /= "NONE" .AND. &
1253 trim(model_key) /= trim(xc_fun_key))
1254 use_skala_model = (index(trim(model_key), "SKALA") > 0)
1255 lwd_kernel_key = adjustl(params%lwd_kernel)
1256 CALL uppercase(lwd_kernel_key)
1257 IF (trim(lwd_kernel_key) == "AUTO") THEN
1258 params%lwd_kernel = "DEFAULT"
1259#if defined(GAUXC_HAS_CUTLASS)
1260 IF (use_skala_model .AND. trim(params%int_exec_space) == "DEVICE") THEN
1261 params%lwd_kernel = "SCHEME1-CUTLASS"
1262 END IF
1263#endif
1264 END IF
1265 params%model_eval_name = model_name
1266 IF (.NOT. params%use_gauxc_model) THEN
1267 ! MODEL NONE and MODEL equal to FUNCTIONAL select conventional GauXC.
1268 params%model_eval_name = "NONE"
1269 END IF
1270 IF (gapw_pseudopotentials .AND. params%use_gauxc_model .AND. .NOT. dft_control%qs_control%gapw_xc .AND. &
1271 .NOT. gapw_paw_pseudopotentials .AND. para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
1272 IF (scf_env%iter_count == 1) THEN
1273 CALL cp_warn( &
1274 __location__, &
1275 "GauXC Skala with METHOD GAPW and GPW_TYPE pseudopotentials evaluates "// &
1276 "the XC term directly on the molecular AO/valence density; no GAPW one-center "// &
1277 "XC correction is used for those regular-grid kinds.")
1278 END IF
1279 END IF
1280 IF (params%device_runtime_fill_fraction <= 0.0_dp .OR. params%device_runtime_fill_fraction > 1.0_dp) THEN
1281 CALL cp_abort(__location__, &
1282 "GAUXC%DEVICE_RUNTIME_FILL_FRACTION must be > 0 and <= 1.")
1283 END IF
1284 IF (atom_chunk_size < -1) THEN
1285 CALL cp_abort(__location__, &
1286 "GAUXC%MODEL_ATOM_CHUNK_SIZE must be -1, zero, or positive.")
1287 END IF
1288 IF (molecular_virial_debug) THEN
1289 IF (molecular_virial_debug_dx <= 0.0_dp) THEN
1290 CALL cp_abort(__location__, &
1291 "GauXC MOLECULAR_VIRIAL_DEBUG_DX must be positive.")
1292 END IF
1293 molecular_virial = .true.
1294 END IF
1295 need_xc_gradient = calculate_forces .OR. molecular_virial
1296 IF (low_spin_scope) THEN
1297 IF (.NOT. dft_control%qs_control%gapw .OR. dft_control%qs_control%gapw_xc) THEN
1298 cpabort("GauXC low spin ROKS requires METHOD GAPW.")
1299 END IF
1300 IF (.NOT. gauxc_gapw_all_all_electron(qs_kind_set)) THEN
1301 cpabort("GAPW low spin ROKS GauXC requires all-electron kinds.")
1302 END IF
1303 IF (is_periodic .OR. do_kpoints .OR. nimages /= 1) THEN
1304 cpabort("GAPW low spin ROKS GauXC requires an isolated Gamma-point calculation.")
1305 END IF
1306 IF (need_xc_gradient) THEN
1307 CALL cp_abort(__location__, &
1308 "GAPW low spin ROKS GauXC nuclear gradients and molecular virial are not implemented.")
1309 END IF
1310 IF (trim(int_exec_space_key) /= "HOST" .OR. trim(params%lb_exec_space) /= "HOST") THEN
1311 cpabort("GAPW low spin ROKS GauXC currently requires HOST execution.")
1312 END IF
1313 IF (params%use_gauxc_model .AND. .NOT. use_skala_model) THEN
1314 cpabort("GAPW low spin ROKS GauXC supports only the SKALA model or conventional XC.")
1315 END IF
1316 IF (use_skala_model .AND. trim(skala_runtime_key) == "MPI") THEN
1317 cpabort("GAPW low spin ROKS Skala requires SKALA_RUNTIME AUTO or SELF.")
1318 END IF
1319 END IF
1320 CALL ensure_gauxc_periodic_reference_scope( &
1321 dft_control, cell, qs_kind_set, do_kpoints, periodic_reference)
1322 IF (is_periodic .AND. periodic_reference .AND. para_env%mepos == 0) THEN
1323 IF (ASSOCIATED(scf_env)) THEN
1324 IF (scf_env%iter_count == 1) THEN
1325 CALL cp_warn( &
1326 __location__, &
1327 "GAUXC%PERIODIC_REFERENCE uses GauXC molecular quadrature for isolated validation "// &
1328 "cells. Compact periodic materials require a dedicated periodic GauXC interface.")
1329 END IF
1330 END IF
1331 END IF
1332 IF (params%use_gauxc_model) THEN
1333 IF (has_nlcc(qs_kind_set)) THEN
1334 CALL cp_abort(__location__, &
1335 "GauXC Skala with NLCC pseudopotentials is not implemented. "// &
1336 "The frozen core density would need a SKALA-consistent feature definition.")
1337 END IF
1338 END IF
1339 IF (params%use_gauxc_model) THEN
1340 CALL set_gauxc_model_atom_chunk_env( &
1341 atom_chunk_size, atom_chunk_size_explicit)
1342 IF (.NOT. grid_explicit) params%grid_type = "SUPERFINE"
1343 IF (.NOT. pruning_explicit) params%pruning_scheme = "UNPRUNED"
1344
1345 grid_key = adjustl(params%grid_type)
1346 pruning_key = adjustl(params%pruning_scheme)
1347 CALL uppercase(grid_key)
1348 CALL uppercase(pruning_key)
1349 IF (use_skala_model .AND. need_xc_gradient .AND. &
1350 (trim(grid_key) /= "SUPERFINE" .OR. trim(pruning_key) /= "UNPRUNED")) THEN
1351 CALL cp_warn( &
1352 __location__, &
1353 "GauXC Skala nuclear gradients are sensitive to the GauXC molecular grid. "// &
1354 "Use GRID SUPERFINE and PRUNING_SCHEME UNPRUNED for quantitative force checks.")
1355 END IF
1356 IF (trim(model_key) == "SKALA") THEN
1357 model_name = ""
1358 env_status = 1
1359 IF (trim(int_exec_space_key) == "DEVICE") THEN
1360 CALL get_environment_variable("SKALA_CUDA_MODEL", model_name, status=env_status)
1361 END IF
1362 IF (env_status /= 0 .OR. len_trim(model_name) == 0) THEN
1363 CALL get_environment_variable("SKALA_MODEL", model_name, status=env_status)
1364 END IF
1365 IF (env_status /= 0 .OR. len_trim(model_name) == 0) THEN
1366 IF (trim(int_exec_space_key) == "DEVICE") THEN
1367 CALL cp_abort( &
1368 __location__, &
1369 "MODEL SKALA with DEVICE execution requires SKALA_CUDA_MODEL or "// &
1370 "GAUXC_SKALA_MODEL")
1371 ELSE
1372 cpabort("MODEL SKALA requires the GAUXC_SKALA_MODEL environment variable")
1373 END IF
1374 END IF
1375 params%model_eval_name = model_name
1376 END IF
1377 END IF
1378 SELECT CASE (trim(skala_runtime_key))
1379 CASE ("AUTO")
1380 params%use_self_runtime = use_skala_model .AND. para_env%num_pe > 1 .AND. params%nspins > 1
1381 CASE ("MPI")
1382 params%use_self_runtime = .false.
1383 CASE ("SELF")
1384 params%use_self_runtime = use_skala_model .AND. para_env%num_pe > 1
1385 CASE DEFAULT
1386 CALL cp_abort(__location__, "Unknown GAUXC%SKALA_RUNTIME value.")
1387 END SELECT
1388 IF (.NOT. use_skala_model) params%use_self_runtime = .false.
1389 params%gradient_runtime = gradient_runtime_key
1390 IF (use_skala_model .AND. para_env%num_pe > 1 .AND. .NOT. params%use_self_runtime .AND. &
1391 para_env%mepos == 0 .AND. ASSOCIATED(scf_env)) THEN
1392 IF (scf_env%iter_count == 1) THEN
1393 CALL cp_warn( &
1394 __location__, &
1395 "GAUXC%SKALA_RUNTIME uses the MPI communicator for energy/VXC. "// &
1396 "SKALA Torch atom chunks can be distributed across MPI ranks; "// &
1397 "set GAUXC_ONEDFT_DISTRIBUTED_TORCH=0 to force rank-0 Torch inference.")
1398 END IF
1399 END IF
1400
1401 ! Resolve the gradient route once the cache knows the actual basis angular momentum.
1402 params%gapw_method = gapw_method
1403 params%need_xc_gradient = need_xc_gradient
1404
1405 IF (.NOT. ASSOCIATED(qs_env%gauxc_cache)) ALLOCATE (qs_env%gauxc_cache)
1406 cache => qs_env%gauxc_cache
1407 CALL gauxc_cache_init( &
1408 cache, &
1409 params, &
1410 para_env, &
1411 particle_set, &
1412 qs_kind_set, &
1413 gauxc_status)
1414
1415 hdf5_output = (trim(output_path) /= "")
1416 write_hdf5_output = hdf5_output .AND. para_env%mepos == 0 .AND. .NOT. external_density
1417 IF (write_hdf5_output .AND. ASSOCIATED(scf_env)) THEN
1418 write_hdf5_output = scf_env%iter_count == 1
1419 END IF
1420 IF (write_hdf5_output) THEN
1421 CALL gauxc_write_molecule_hdf5( &
1422 cache%molecule, &
1423 output_path, &
1424 "molecule.h5", &
1425 "molecule", &
1426 gauxc_status)
1427 CALL gauxc_check_status(gauxc_status)
1428 CALL gauxc_write_basisset_hdf5( &
1429 cache%basisset, &
1430 output_path, &
1431 "basisset.h5", &
1432 "basisset", &
1433 gauxc_status)
1434 CALL gauxc_check_status(gauxc_status)
1435 END IF
1436
1437 IF (qs_env%run_rtp) THEN
1438 cpabort("GAUXC XC energy currently does not support real-time propagation")
1439 END IF
1440
1441 exc_value = 0.0_dp
1442
1443 IF (ASSOCIATED(matrix_vxc)) CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1444 CALL dbcsr_allocate_matrix_set(matrix_vxc, params%nspins)
1445
1446 DO img = 1, nimages
1447 IF (img > 1) THEN
1448 cpabort("UNIMPLEMENTED: Handling nimg>1 in k-point integration")
1449 END IF
1450 CALL dbcsr_to_dense(rho_ao(1, img), density_scalar, para_env)
1451 CALL para_env%sum(density_scalar)
1452 IF (params%nspins == 1) THEN
1453 gauxc_xc_result = gauxc_compute_xc( &
1454 cache%integrator, &
1455 density_scalar, &
1456 nspins=params%nspins, &
1457 status=gauxc_status, &
1458 model=trim(params%model_eval_name))
1459 CALL gauxc_check_status(gauxc_status)
1460 IF (need_xc_gradient) THEN
1461 IF (params%use_fd_gradient) THEN
1462 CALL gauxc_xc_gradient_fd( &
1463 particle_set, qs_kind_set, density_scalar, params%nspins, params%model_eval_name, &
1464 params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1465 params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1466 params%device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
1467 exc_grad%exc_grad)
1468 ELSE IF (params%use_gradient_self_runtime) THEN
1469 exc_grad = gauxc_compute_xc_gradient( &
1470 cache%gradient_integrator, &
1471 density_scalar, &
1472 nspins=params%nspins, &
1473 natom=params%natom, &
1474 status=gauxc_status, &
1475 model=trim(params%model_eval_name))
1476 ELSE
1477 exc_grad = gauxc_compute_xc_gradient( &
1478 cache%integrator, &
1479 density_scalar, &
1480 nspins=params%nspins, &
1481 natom=params%natom, &
1482 status=gauxc_status, &
1483 model=trim(params%model_eval_name))
1484 END IF
1485 CALL gauxc_check_status(gauxc_status)
1486 IF (calculate_forces) THEN
1487 CALL add_gauxc_gradient_to_force( &
1488 exc_grad%exc_grad, &
1489 force, &
1490 atomic_kind_set, &
1491 para_env)
1492 END IF
1493 IF (molecular_virial) THEN
1494 CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1495 END IF
1496 IF (molecular_virial_debug) THEN
1497 CALL debug_gauxc_molecular_virial( &
1498 exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, params%nspins, &
1499 params%model_eval_name, params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1500 params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1501 params%device_runtime_fill_fraction, molecular_virial_debug_dx, para_env)
1502 END IF
1503 DEALLOCATE (exc_grad%exc_grad)
1504 END IF
1505 ELSE
1506 cpassert(params%nspins == 2)
1507 ! In here:
1508 ! scalar <- rho_ao(1, :) + rho_ao(2, :)
1509 ! zeta <- rho_ao(1, :) - rho_ao(2, :)
1510 CALL dbcsr_to_dense(rho_ao(2, img), density_zeta, para_env)
1511 CALL para_env%sum(density_zeta)
1512 ! Do NOT reorder the following lines!
1513 density_scalar(:, :) = density_scalar(:, :) + density_zeta(:, :)
1514 ! Factor two because the next line is evaluated after the above line.
1515 ! We need to subtract density_zeta once to undo the above line and
1516 ! a second time because that is what UKS requires.
1517 ! This style lowers memory footprint.
1518 density_zeta(:, :) = density_scalar(:, :) - 2.0_dp*density_zeta(:, :)
1519 gauxc_xc_result = gauxc_compute_xc( &
1520 cache%integrator, &
1521 density_scalar, &
1522 density_zeta, &
1523 params%nspins, &
1524 gauxc_status, &
1525 model=trim(params%model_eval_name))
1526 CALL gauxc_check_status(gauxc_status)
1527 IF (need_xc_gradient) THEN
1528 IF (params%use_fd_gradient) THEN
1529 CALL gauxc_xc_gradient_fd( &
1530 particle_set, qs_kind_set, density_scalar, params%nspins, params%model_eval_name, &
1531 params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1532 params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1533 params%device_runtime_fill_fraction, gapw_fd_gradient_dx, para_env, &
1534 exc_grad%exc_grad, density_zeta=density_zeta)
1535 ELSE IF (params%use_gradient_self_runtime) THEN
1536 exc_grad = gauxc_compute_xc_gradient( &
1537 cache%gradient_integrator, &
1538 density_scalar, &
1539 density_zeta, &
1540 params%nspins, &
1541 params%natom, &
1542 gauxc_status, &
1543 model=trim(params%model_eval_name))
1544 ELSE
1545 exc_grad = gauxc_compute_xc_gradient( &
1546 cache%integrator, &
1547 density_scalar, &
1548 density_zeta, &
1549 params%nspins, &
1550 params%natom, &
1551 gauxc_status, &
1552 model=trim(params%model_eval_name))
1553 END IF
1554 CALL gauxc_check_status(gauxc_status)
1555 IF (calculate_forces) THEN
1556 CALL add_gauxc_gradient_to_force( &
1557 exc_grad%exc_grad, &
1558 force, &
1559 atomic_kind_set, &
1560 para_env)
1561 END IF
1562 IF (molecular_virial) THEN
1563 CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1564 END IF
1565 IF (molecular_virial_debug) THEN
1566 CALL debug_gauxc_molecular_virial( &
1567 exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, params%nspins, &
1568 params%model_eval_name, params%xc_fun_name, params%grid_type, params%radial_quadrature, params%pruning_scheme, &
1569 params%lb_exec_space, params%int_exec_space, params%lwd_kernel, params%batch_size, &
1570 params%device_runtime_fill_fraction, molecular_virial_debug_dx, para_env, &
1571 density_zeta=density_zeta)
1572 END IF
1573 DEALLOCATE (exc_grad%exc_grad)
1574 END IF
1575 END IF
1576
1577 exc_value = exc_value + gauxc_xc_result%exc
1578
1579 IF (params%nspins == 1) THEN
1580 IF (img == 1) THEN
1581 matrix_vxc(1) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(1, img))
1582 ELSE
1583 cpabort("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1584 END IF
1585 ELSE
1586 cpassert(params%nspins == 2)
1587 ! Transform derivatives from total/spin density back to alpha/beta channels.
1588 vxc_zeta_tmp = dense_to_dbcsr(gauxc_xc_result%vxc_zeta, rho_ao(1, img))
1589 IF (img == 1) THEN
1590 DO ispin = 1, 2
1591 matrix_vxc(ispin) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(ispin, 1))
1592 CALL dbcsr_add( &
1593 matrix_vxc(ispin)%matrix, &
1594 vxc_zeta_tmp%matrix, &
1595 1.0_dp, &
1596 ! 1.0 for ispin==1, -1.0 for ispin==2
1597 1.0_dp - real(ispin - 1, dp)*2.0_dp)
1598 END DO
1599 ELSE
1600 cpabort("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1601 END IF
1602 CALL dbcsr_release(vxc_zeta_tmp%matrix)
1603 DEALLOCATE (vxc_zeta_tmp%matrix)
1604 END IF
1605 END DO
1606
1607 DEALLOCATE (density_scalar)
1608 IF (ALLOCATED(density_zeta)) DEALLOCATE (density_zeta)
1609 DEALLOCATE (gauxc_xc_result%vxc_scalar)
1610 IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
1611
1612 IF (external_density) THEN
1613 exc_external = exc_value
1614 matrix_vxc_external => matrix_vxc
1615 ELSE
1616 energy%exc = exc_value
1617 CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
1618 END IF
1619 DO ispin = 1, params%nspins
1620 CALL dbcsr_finalize(matrix_vxc(ispin)%matrix)
1621 END DO
1622
1623 END SUBROUTINE apply_gauxc
1624
1625END 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:95
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition cp_dbcsr_operations.F:103
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:637
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:817
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:545
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:4291
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:122
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, rho_ao_external, exc_external, matrix_vxc_external)
Evaluate molecular GauXC, optionally for an explicit configuration density.
Definition xc_gauxc_functional.F:1033
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:787
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:75
external_potential_types::gth_potential_type
Definition external_potential_types.F:121
external_potential_types::sgp_potential_type
Definition external_potential_types.F:179
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:768
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:224
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:100
xc_gauxc_cache::cp_gauxc_cache_params
Definition xc_gauxc_cache.F:33
xc_gauxc_cache::cp_gauxc_cache_type
Definition xc_gauxc_cache.F:58
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