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