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