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