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