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 cp_control_types, ONLY: dft_control_type
12 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
13 : dbcsr_create,&
14 : dbcsr_finalize,&
15 : dbcsr_get_info,&
16 : dbcsr_p_type,&
17 : dbcsr_release
18 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
19 : dbcsr_deallocate_matrix_set
20 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
21 : USE external_potential_types, ONLY: gth_potential_type,&
22 : sgp_potential_type
23 : USE input_constants, ONLY: xc_vdw_fun_nonloc
24 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
25 : section_vals_get_subs_vals2,&
26 : section_vals_type,&
27 : section_vals_val_get
28 : USE iso_c_binding, ONLY: c_double
29 : USE kinds, ONLY: default_path_length,&
30 : default_string_length,&
31 : dp
32 : USE message_passing, ONLY: mp_comm_self,&
33 : mp_para_env_type
34 : USE particle_types, ONLY: particle_type
35 : USE qs_energy_types, ONLY: qs_energy_type
36 : USE qs_environment_types, ONLY: get_qs_env,&
37 : qs_environment_type
38 : USE qs_force_types, ONLY: qs_force_type
39 : USE qs_kind_types, ONLY: get_qs_kind,&
40 : qs_kind_type
41 : USE qs_ks_types, ONLY: qs_ks_env_type,&
42 : set_ks_env
43 : USE qs_rho_types, ONLY: qs_rho_get,&
44 : qs_rho_type
45 : USE qs_scf_types, ONLY: qs_scf_env_type
46 : USE string_utilities, ONLY: uppercase
47 : USE xc_gauxc_interface, ONLY: &
48 : cp_gauxc_basisset_type, cp_gauxc_grid_type, cp_gauxc_integrator_type, &
49 : cp_gauxc_molecule_type, cp_gauxc_status_type, cp_gauxc_xc_gradient_type, cp_gauxc_xc_type, &
50 : gauxc_check_status, gauxc_compute_xc, gauxc_compute_xc_gradient, gauxc_create_basisset, &
51 : gauxc_create_grid, gauxc_create_integrator, gauxc_create_molecule, gauxc_destroy_basisset, &
52 : gauxc_destroy_grid, gauxc_destroy_integrator, gauxc_destroy_molecule, &
53 : gauxc_write_basisset_hdf5, gauxc_write_molecule_hdf5
54 : USE xc_rho_cflags_types, ONLY: xc_rho_cflags_type
55 : #include "../base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 :
61 : LOGICAL, PARAMETER :: debug_this_module = .TRUE.
62 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_gauxc_functional'
63 :
64 : PUBLIC :: apply_gauxc, skala_info, xc_section_uses_gauxc
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief ...
70 : !> \param dbcsr_mat ...
71 : !> \param dense_mat ...
72 : ! **************************************************************************************************
73 636 : SUBROUTINE dbcsr_to_dense(dbcsr_mat, dense_mat)
74 : USE cp_dbcsr_api, ONLY: dbcsr_get_info, dbcsr_get_matrix_type, dbcsr_iterator_type, &
75 : dbcsr_iterator_start, dbcsr_iterator_next_block, dbcsr_iterator_blocks_left, &
76 : dbcsr_iterator_stop, dbcsr_type_antisymmetric, dbcsr_type_symmetric
77 : TYPE(dbcsr_p_type), INTENT(IN) :: dbcsr_mat
78 : REAL(c_double), ALLOCATABLE, DIMENSION(:, :), &
79 : INTENT(INOUT) :: dense_mat
80 :
81 : CHARACTER :: matrix_type
82 : INTEGER :: col, col_end, col_start, icol, irow, &
83 : nblkcols_total, nblkrows_total, ncols, &
84 : nrows, row, row_end, row_start
85 636 : INTEGER, ALLOCATABLE, DIMENSION(:) :: c_offset, r_offset
86 636 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
87 : LOGICAL :: transposed
88 636 : REAL(c_double), POINTER :: block(:, :)
89 : TYPE(dbcsr_iterator_type) :: iter
90 :
91 : CALL dbcsr_get_info(dbcsr_mat%matrix, &
92 : row_blk_size=row_blk_size, &
93 : col_blk_size=col_blk_size, &
94 : nblkrows_total=nblkrows_total, &
95 : nblkcols_total=nblkcols_total, &
96 : nfullrows_total=nrows, &
97 636 : nfullcols_total=ncols)
98 636 : matrix_type = dbcsr_get_matrix_type(dbcsr_mat%matrix)
99 :
100 636 : IF (.NOT. ALLOCATED(dense_mat)) THEN
101 2544 : ALLOCATE (dense_mat(nrows, ncols))
102 0 : ELSE IF (.NOT. ALL(SHAPE(dense_mat) == [nrows, ncols])) THEN
103 0 : DEALLOCATE (dense_mat)
104 0 : ALLOCATE (dense_mat(nrows, ncols))
105 : ELSE
106 0 : CPASSERT(ALL(SHAPE(dense_mat) == [nrows, ncols]))
107 : END IF
108 636 : dense_mat = 0._dp
109 :
110 3180 : ALLOCATE (r_offset(nblkrows_total), c_offset(nblkcols_total))
111 :
112 636 : r_offset(1) = 1
113 1448 : DO row = 2, nblkrows_total
114 1448 : r_offset(row) = r_offset(row - 1) + row_blk_size(row - 1)
115 : END DO
116 636 : c_offset(1) = 1
117 1448 : DO col = 2, nblkcols_total
118 1448 : c_offset(col) = c_offset(col - 1) + col_blk_size(col - 1)
119 : END DO
120 :
121 636 : CALL dbcsr_iterator_start(iter, dbcsr_mat%matrix)
122 2060 : DO WHILE (dbcsr_iterator_blocks_left(iter))
123 1424 : CALL dbcsr_iterator_next_block(iter, irow, icol, block, transposed=transposed)
124 1424 : row_start = r_offset(irow)
125 1424 : row_end = row_start + row_blk_size(irow) - 1
126 1424 : col_start = c_offset(icol)
127 1424 : col_end = col_start + col_blk_size(icol) - 1
128 2060 : IF (transposed) THEN
129 0 : dense_mat(row_start:row_end, col_start:col_end) = TRANSPOSE(block)
130 0 : IF (irow /= icol) THEN
131 0 : IF (matrix_type == dbcsr_type_symmetric) THEN
132 0 : dense_mat(col_start:col_end, row_start:row_end) = block
133 0 : ELSE IF (matrix_type == dbcsr_type_antisymmetric) THEN
134 0 : dense_mat(col_start:col_end, row_start:row_end) = -block
135 : END IF
136 : END IF
137 : ELSE
138 71224 : dense_mat(row_start:row_end, col_start:col_end) = block
139 1424 : IF (irow /= icol) THEN
140 628 : IF (matrix_type == dbcsr_type_symmetric) THEN
141 25332 : dense_mat(col_start:col_end, row_start:row_end) = TRANSPOSE(block)
142 0 : ELSE IF (matrix_type == dbcsr_type_antisymmetric) THEN
143 0 : dense_mat(col_start:col_end, row_start:row_end) = -TRANSPOSE(block)
144 : END IF
145 : END IF
146 : END IF
147 : END DO
148 636 : CALL dbcsr_iterator_stop(iter)
149 :
150 636 : DEALLOCATE (r_offset, c_offset)
151 :
152 1272 : END SUBROUTINE dbcsr_to_dense
153 :
154 : ! ******, ***********************************************************************************
155 : !> \brief Convert a dense symmetric matrix to a DBCSR matrix with full upper block structure.
156 : !> This creates all upper-triangular blocks, not just those present in a template.
157 : !> This is needed because GauXC computes VXC for the full dense density matrix.
158 : !> \param dense_mat Input dense matrix
159 : !> \param template_dbcsr Template DBCSR matrix for distribution and block sizes
160 : !> \return dbcsr_mat Output DBCSR matrix with full upper block structure
161 : ! **************************************************************************************************
162 1524 : FUNCTION dense_to_dbcsr(dense_mat, template_dbcsr) RESULT(dbcsr_mat)
163 : USE cp_dbcsr_api, ONLY: &
164 : dbcsr_create, &
165 : dbcsr_distribution_get, &
166 : dbcsr_distribution_type, &
167 : dbcsr_finalize, &
168 : dbcsr_get_info, &
169 : dbcsr_get_stored_coordinates, &
170 : dbcsr_init_p, &
171 : dbcsr_put_block, &
172 : dbcsr_release, &
173 : dbcsr_type_symmetric, &
174 : dbcsr_work_create
175 : REAL(c_double), DIMENSION(:, :), INTENT(IN) :: dense_mat
176 : TYPE(dbcsr_p_type), INTENT(IN) :: template_dbcsr
177 : TYPE(dbcsr_p_type) :: dbcsr_mat
178 :
179 : INTEGER :: col, icol, irow, mynode, nblkcols_total, &
180 : nblkrows_total, ncols, nrows, owner, &
181 : row
182 762 : INTEGER, ALLOCATABLE, DIMENSION(:) :: c_offset, r_offset
183 762 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
184 : TYPE(dbcsr_distribution_type) :: dist
185 :
186 : CALL dbcsr_get_info(template_dbcsr%matrix, &
187 : row_blk_size=row_blk_size, &
188 : col_blk_size=col_blk_size, &
189 : nblkrows_total=nblkrows_total, &
190 : nblkcols_total=nblkcols_total, &
191 : nfullrows_total=nrows, &
192 : nfullcols_total=ncols, &
193 762 : distribution=dist)
194 762 : CALL dbcsr_distribution_get(dist, mynode=mynode)
195 :
196 762 : CPASSERT(nrows == SIZE(dense_mat, 1))
197 762 : CPASSERT(ncols == SIZE(dense_mat, 2))
198 :
199 762 : CALL dbcsr_init_p(dbcsr_mat%matrix)
200 : CALL dbcsr_create(dbcsr_mat%matrix, &
201 : template=template_dbcsr%matrix, &
202 : name="VXC from GauXC (dense)", &
203 762 : matrix_type=dbcsr_type_symmetric)
204 762 : CALL dbcsr_work_create(dbcsr_mat%matrix, work_mutable=.TRUE.)
205 :
206 3810 : ALLOCATE (r_offset(nblkrows_total), c_offset(nblkcols_total))
207 :
208 762 : r_offset(1) = 1
209 1680 : DO row = 2, nblkrows_total
210 1680 : r_offset(row) = r_offset(row - 1) + row_blk_size(row - 1)
211 : END DO
212 762 : c_offset(1) = 1
213 1680 : DO col = 2, nblkcols_total
214 1680 : c_offset(col) = c_offset(col - 1) + col_blk_size(col - 1)
215 : END DO
216 :
217 2442 : DO irow = 1, nblkrows_total
218 6702 : DO icol = 1, nblkcols_total
219 4260 : IF (irow > icol) CYCLE
220 2970 : CALL dbcsr_get_stored_coordinates(dbcsr_mat%matrix, irow, icol, owner)
221 2970 : IF (owner /= mynode) CYCLE
222 : CALL dbcsr_put_block(dbcsr_mat%matrix, irow, icol, &
223 : 0.5_dp*( &
224 : dense_mat(r_offset(irow):r_offset(irow) + row_blk_size(irow) - 1, &
225 : c_offset(icol):c_offset(icol) + col_blk_size(icol) - 1) + &
226 : TRANSPOSE(dense_mat(r_offset(icol):r_offset(icol) + row_blk_size(icol) - 1, &
227 89795 : c_offset(irow):c_offset(irow) + col_blk_size(irow) - 1))))
228 : END DO
229 : END DO
230 :
231 762 : CALL dbcsr_finalize(dbcsr_mat%matrix)
232 :
233 762 : DEALLOCATE (r_offset, c_offset)
234 :
235 762 : END FUNCTION dense_to_dbcsr
236 :
237 : ! **************************************************************************************************
238 : !> \brief ...
239 : !> \param xc_section ...
240 : !> \return ...
241 : ! **************************************************************************************************
242 510 : FUNCTION get_gauxc_functional(xc_section) RESULT(gauxc_functional_section)
243 : TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
244 : TYPE(section_vals_type), POINTER :: gauxc_functional_section
245 :
246 : INTEGER :: ifun
247 : TYPE(section_vals_type), POINTER :: functionals, xc_fun
248 :
249 510 : NULLIFY (gauxc_functional_section)
250 :
251 510 : functionals => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
252 510 : IF (.NOT. ASSOCIATED(functionals)) THEN
253 0 : CPABORT("XC_FUNCTIONAL section not found")
254 : END IF
255 :
256 510 : ifun = 0
257 : DO
258 1020 : ifun = ifun + 1
259 1020 : xc_fun => section_vals_get_subs_vals2(functionals, i_section=ifun)
260 1020 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
261 510 : IF (xc_fun%section%name /= "GAUXC" .OR. ifun > 1) THEN
262 0 : CPABORT("GauXC functionals are mutually exclusive with any other functional.")
263 : END IF
264 510 : gauxc_functional_section => xc_fun
265 : END DO
266 :
267 510 : IF (.NOT. ASSOCIATED(gauxc_functional_section)) THEN
268 0 : CPABORT("No XC functional found in XC_FUNCTIONAL section")
269 : END IF
270 510 : END FUNCTION get_gauxc_functional
271 :
272 : ! **************************************************************************************************
273 : !> \brief ...
274 : !> \param xc_section ...
275 : !> \return ...
276 : ! **************************************************************************************************
277 14128 : FUNCTION xc_section_uses_gauxc(xc_section) RESULT(uses_gauxc)
278 : TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
279 : LOGICAL :: uses_gauxc
280 :
281 : INTEGER :: ifun
282 : TYPE(section_vals_type), POINTER :: functionals, xc_fun
283 :
284 14128 : uses_gauxc = .FALSE.
285 14128 : IF (.NOT. ASSOCIATED(xc_section)) RETURN
286 :
287 14128 : functionals => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
288 14128 : IF (.NOT. ASSOCIATED(functionals)) RETURN
289 :
290 14128 : ifun = 0
291 : DO
292 27598 : ifun = ifun + 1
293 27598 : xc_fun => section_vals_get_subs_vals2(functionals, i_section=ifun)
294 27598 : IF (.NOT. ASSOCIATED(xc_fun)) EXIT
295 27598 : IF (xc_fun%section%name == "GAUXC") THEN
296 : uses_gauxc = .TRUE.
297 : EXIT
298 : END IF
299 : END DO
300 :
301 : END FUNCTION xc_section_uses_gauxc
302 :
303 : ! **************************************************************************************************
304 : !> \brief Reject unsupported pseudopotential variants in GauXC GAPW mode
305 : !> \param qs_kind_set ...
306 : ! **************************************************************************************************
307 52 : SUBROUTINE ensure_gauxc_gapw_all_electron(qs_kind_set)
308 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
309 :
310 : INTEGER :: ikind
311 : TYPE(gth_potential_type), POINTER :: gth_potential
312 : TYPE(sgp_potential_type), POINTER :: sgp_potential
313 :
314 52 : CPASSERT(ASSOCIATED(qs_kind_set))
315 :
316 104 : DO ikind = 1, SIZE(qs_kind_set)
317 52 : NULLIFY (gth_potential, sgp_potential)
318 : CALL get_qs_kind(qs_kind_set(ikind), &
319 : gth_potential=gth_potential, &
320 52 : sgp_potential=sgp_potential)
321 104 : IF (ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)) THEN
322 : CALL cp_abort(__LOCATION__, &
323 : "GauXC with METHOD GAPW currently supports all-electron potentials only. "// &
324 0 : "Use POTENTIAL ALL for GAPW validation or METHOD GPW with pseudopotentials.")
325 : END IF
326 : END DO
327 :
328 52 : END SUBROUTINE ensure_gauxc_gapw_all_electron
329 :
330 : ! **************************************************************************************************
331 : !> \brief adds a replicated GauXC energy gradient to the local CP2K force accumulator
332 : !> \param exc_grad ...
333 : !> \param force ...
334 : !> \param atomic_kind_set ...
335 : !> \param para_env ...
336 : ! **************************************************************************************************
337 4 : SUBROUTINE add_gauxc_gradient_to_force(exc_grad, force, atomic_kind_set, para_env)
338 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
339 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
340 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
341 : TYPE(mp_para_env_type), POINTER :: para_env
342 :
343 : INTEGER :: ia, iatom, ikind, natom_kind
344 : TYPE(atomic_kind_type), POINTER :: atomic_kind
345 :
346 4 : CPASSERT(ASSOCIATED(force))
347 4 : CPASSERT(ASSOCIATED(atomic_kind_set))
348 :
349 4 : IF (para_env%mepos /= 0) RETURN
350 :
351 5 : DO ikind = 1, SIZE(atomic_kind_set, 1)
352 3 : atomic_kind => atomic_kind_set(ikind)
353 3 : CALL get_atomic_kind(atomic_kind=atomic_kind, natom=natom_kind)
354 11 : DO ia = 1, natom_kind
355 6 : iatom = atomic_kind%atom_list(ia)
356 : force(ikind)%rho_elec(:, ia) = force(ikind)%rho_elec(:, ia) + &
357 27 : exc_grad(3*iatom - 2:3*iatom)
358 : END DO
359 : END DO
360 :
361 : END SUBROUTINE add_gauxc_gradient_to_force
362 :
363 : ! **************************************************************************************************
364 : !> \brief compute a GauXC XC energy for diagnostic finite differences
365 : !> \param particle_set_eval ...
366 : !> \param qs_kind_set ...
367 : !> \param density_scalar ...
368 : !> \param nspins ...
369 : !> \param model_name ...
370 : !> \param xc_fun_name ...
371 : !> \param grid_type ...
372 : !> \param radial_quadrature ...
373 : !> \param pruning_scheme ...
374 : !> \param lb_exec_space ...
375 : !> \param int_exec_space ...
376 : !> \param batch_size ...
377 : !> \param exc ...
378 : !> \param density_zeta ...
379 : ! **************************************************************************************************
380 0 : SUBROUTINE gauxc_xc_energy_for_particles( &
381 0 : particle_set_eval, qs_kind_set, density_scalar, nspins, model_name, &
382 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
383 0 : int_exec_space, batch_size, exc, density_zeta)
384 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set_eval
385 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
386 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
387 : INTEGER, INTENT(IN) :: nspins
388 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, &
389 : radial_quadrature, pruning_scheme, &
390 : lb_exec_space, int_exec_space
391 : INTEGER, INTENT(IN) :: batch_size
392 : REAL(KIND=dp), INTENT(OUT) :: exc
393 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
394 : OPTIONAL :: density_zeta
395 :
396 : TYPE(cp_gauxc_basisset_type) :: gauxc_basis_fd
397 : TYPE(cp_gauxc_grid_type) :: gauxc_grid_fd
398 : TYPE(cp_gauxc_integrator_type) :: gauxc_integrator_fd
399 : TYPE(cp_gauxc_molecule_type) :: gauxc_mol_fd
400 : TYPE(cp_gauxc_status_type) :: gauxc_status
401 0 : TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
402 :
403 0 : gauxc_mol_fd = gauxc_create_molecule(particle_set_eval, gauxc_status)
404 0 : CALL gauxc_check_status(gauxc_status)
405 0 : gauxc_basis_fd = gauxc_create_basisset(qs_kind_set, particle_set_eval, gauxc_status)
406 0 : CALL gauxc_check_status(gauxc_status)
407 : gauxc_grid_fd = gauxc_create_grid( &
408 : gauxc_mol_fd, &
409 : gauxc_basis_fd, &
410 : grid_type, &
411 : radial_quadrature, &
412 : pruning_scheme, &
413 : lb_exec_space, &
414 : batch_size, &
415 : gauxc_status, &
416 0 : mpi_comm=mp_comm_self%get_handle())
417 0 : CALL gauxc_check_status(gauxc_status)
418 : gauxc_integrator_fd = gauxc_create_integrator( &
419 : TRIM(xc_fun_name), &
420 : gauxc_grid_fd, &
421 : int_exec_space, &
422 : nspins, &
423 0 : gauxc_status)
424 0 : CALL gauxc_check_status(gauxc_status)
425 :
426 0 : IF (nspins == 1) THEN
427 : gauxc_xc_result = gauxc_compute_xc( &
428 : gauxc_integrator_fd, &
429 : density_scalar, &
430 : nspins=nspins, &
431 : status=gauxc_status, &
432 0 : model=TRIM(model_name))
433 : ELSE
434 0 : CPASSERT(nspins == 2)
435 0 : CPASSERT(PRESENT(density_zeta))
436 : gauxc_xc_result = gauxc_compute_xc( &
437 : gauxc_integrator_fd, &
438 : density_scalar, &
439 : density_zeta, &
440 : nspins, &
441 : gauxc_status, &
442 0 : model=TRIM(model_name))
443 : END IF
444 0 : CALL gauxc_check_status(gauxc_status)
445 0 : exc = gauxc_xc_result%exc
446 :
447 0 : IF (ALLOCATED(gauxc_xc_result%vxc_scalar)) DEALLOCATE (gauxc_xc_result%vxc_scalar)
448 0 : IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
449 :
450 0 : CALL gauxc_destroy_integrator(gauxc_integrator_fd, gauxc_status)
451 0 : CALL gauxc_check_status(gauxc_status)
452 0 : CALL gauxc_destroy_grid(gauxc_grid_fd, gauxc_status)
453 0 : CALL gauxc_check_status(gauxc_status)
454 0 : CALL gauxc_destroy_basisset(gauxc_basis_fd, gauxc_status)
455 0 : CALL gauxc_check_status(gauxc_status)
456 0 : CALL gauxc_destroy_molecule(gauxc_mol_fd, gauxc_status)
457 0 : CALL gauxc_check_status(gauxc_status)
458 :
459 0 : END SUBROUTINE gauxc_xc_energy_for_particles
460 :
461 : ! **************************************************************************************************
462 : !> \brief compute a finite-difference GauXC XC nuclear gradient at fixed density
463 : !> \param particle_set ...
464 : !> \param qs_kind_set ...
465 : !> \param density_scalar ...
466 : !> \param nspins ...
467 : !> \param model_name ...
468 : !> \param xc_fun_name ...
469 : !> \param grid_type ...
470 : !> \param radial_quadrature ...
471 : !> \param pruning_scheme ...
472 : !> \param lb_exec_space ...
473 : !> \param int_exec_space ...
474 : !> \param batch_size ...
475 : !> \param dx ...
476 : !> \param para_env ...
477 : !> \param exc_grad ...
478 : !> \param density_zeta ...
479 : ! **************************************************************************************************
480 0 : SUBROUTINE gauxc_xc_gradient_fd( &
481 0 : particle_set, qs_kind_set, density_scalar, nspins, model_name, &
482 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
483 0 : int_exec_space, batch_size, dx, para_env, exc_grad, density_zeta)
484 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
485 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
486 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
487 : INTEGER, INTENT(IN) :: nspins
488 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, &
489 : radial_quadrature, pruning_scheme, &
490 : lb_exec_space, int_exec_space
491 : INTEGER, INTENT(IN) :: batch_size
492 : REAL(KIND=dp), INTENT(IN) :: dx
493 : TYPE(mp_para_env_type), POINTER :: para_env
494 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:), &
495 : INTENT(OUT) :: exc_grad
496 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
497 : OPTIONAL :: density_zeta
498 :
499 : INTEGER :: iatom, idir
500 : REAL(KIND=dp) :: xc_minus, xc_plus
501 0 : TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus
502 :
503 0 : CPASSERT(ASSOCIATED(particle_set))
504 0 : CPASSERT(dx > 0.0_dp)
505 :
506 0 : ALLOCATE (exc_grad(3*SIZE(particle_set)))
507 0 : exc_grad = 0.0_dp
508 :
509 0 : IF (para_env%mepos == 0) THEN
510 0 : ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set)))
511 :
512 0 : DO iatom = 1, SIZE(particle_set)
513 0 : DO idir = 1, 3
514 0 : particle_set_minus = particle_set
515 0 : particle_set_plus = particle_set
516 0 : particle_set_minus(iatom)%r(idir) = particle_set_minus(iatom)%r(idir) - dx
517 0 : particle_set_plus(iatom)%r(idir) = particle_set_plus(iatom)%r(idir) + dx
518 0 : IF (PRESENT(density_zeta)) THEN
519 : CALL gauxc_xc_energy_for_particles( &
520 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
521 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
522 0 : int_exec_space, batch_size, xc_plus, density_zeta=density_zeta)
523 : CALL gauxc_xc_energy_for_particles( &
524 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
525 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
526 0 : int_exec_space, batch_size, xc_minus, density_zeta=density_zeta)
527 : ELSE
528 : CALL gauxc_xc_energy_for_particles( &
529 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
530 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
531 0 : int_exec_space, batch_size, xc_plus)
532 : CALL gauxc_xc_energy_for_particles( &
533 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
534 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
535 0 : int_exec_space, batch_size, xc_minus)
536 : END IF
537 0 : exc_grad(3*iatom - 3 + idir) = (xc_plus - xc_minus)/(2.0_dp*dx)
538 : END DO
539 : END DO
540 :
541 0 : DEALLOCATE (particle_set_minus, particle_set_plus)
542 : END IF
543 :
544 0 : CALL para_env%bcast(exc_grad, 0)
545 :
546 0 : END SUBROUTINE gauxc_xc_gradient_fd
547 :
548 : ! **************************************************************************************************
549 : !> \brief finite-difference check of the molecular GauXC XC virial diagnostic
550 : !> \param exc_grad ...
551 : !> \param particle_set ...
552 : !> \param qs_kind_set ...
553 : !> \param density_scalar ...
554 : !> \param nspins ...
555 : !> \param model_name ...
556 : !> \param xc_fun_name ...
557 : !> \param grid_type ...
558 : !> \param radial_quadrature ...
559 : !> \param pruning_scheme ...
560 : !> \param lb_exec_space ...
561 : !> \param int_exec_space ...
562 : !> \param batch_size ...
563 : !> \param dx ...
564 : !> \param para_env ...
565 : !> \param density_zeta ...
566 : ! **************************************************************************************************
567 0 : SUBROUTINE debug_gauxc_molecular_virial( &
568 0 : exc_grad, particle_set, qs_kind_set, density_scalar, nspins, model_name, &
569 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
570 0 : int_exec_space, batch_size, dx, para_env, density_zeta)
571 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
572 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
573 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
574 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: density_scalar
575 : INTEGER, INTENT(IN) :: nspins
576 : CHARACTER(len=*), INTENT(IN) :: model_name, xc_fun_name, grid_type, &
577 : radial_quadrature, pruning_scheme, &
578 : lb_exec_space, int_exec_space
579 : INTEGER, INTENT(IN) :: batch_size
580 : REAL(KIND=dp), INTENT(IN) :: dx
581 : TYPE(mp_para_env_type), POINTER :: para_env
582 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
583 : OPTIONAL :: density_zeta
584 :
585 : INTEGER :: iatom, iw
586 : REAL(KIND=dp) :: analytic_trace, diff_trace, &
587 : numerical_trace, xc_minus, xc_plus
588 : REAL(KIND=dp), DIMENSION(3) :: center, displacement, grad
589 0 : TYPE(particle_type), ALLOCATABLE, DIMENSION(:) :: particle_set_minus, particle_set_plus
590 :
591 0 : CPASSERT(ASSOCIATED(particle_set))
592 0 : CPASSERT(SIZE(exc_grad) == 3*SIZE(particle_set))
593 :
594 0 : IF (para_env%mepos /= 0) RETURN
595 :
596 0 : center = 0.0_dp
597 0 : DO iatom = 1, SIZE(particle_set)
598 0 : center = center + particle_set(iatom)%r
599 : END DO
600 0 : center = center/REAL(SIZE(particle_set), dp)
601 :
602 0 : ALLOCATE (particle_set_minus(SIZE(particle_set)), particle_set_plus(SIZE(particle_set)))
603 0 : particle_set_minus = particle_set
604 0 : particle_set_plus = particle_set
605 :
606 0 : analytic_trace = 0.0_dp
607 0 : DO iatom = 1, SIZE(particle_set)
608 0 : grad = exc_grad(3*iatom - 2:3*iatom)
609 0 : displacement = particle_set(iatom)%r - center
610 0 : analytic_trace = analytic_trace + DOT_PRODUCT(grad, displacement)
611 0 : particle_set_minus(iatom)%r = center + (1.0_dp - dx)*displacement
612 0 : particle_set_plus(iatom)%r = center + (1.0_dp + dx)*displacement
613 : END DO
614 0 : analytic_trace = analytic_trace/3.0_dp
615 :
616 0 : IF (PRESENT(density_zeta)) THEN
617 : CALL gauxc_xc_energy_for_particles( &
618 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
619 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
620 0 : int_exec_space, batch_size, xc_plus, density_zeta=density_zeta)
621 : CALL gauxc_xc_energy_for_particles( &
622 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
623 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
624 0 : int_exec_space, batch_size, xc_minus, density_zeta=density_zeta)
625 : ELSE
626 : CALL gauxc_xc_energy_for_particles( &
627 : particle_set_plus, qs_kind_set, density_scalar, nspins, model_name, &
628 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
629 0 : int_exec_space, batch_size, xc_plus)
630 : CALL gauxc_xc_energy_for_particles( &
631 : particle_set_minus, qs_kind_set, density_scalar, nspins, model_name, &
632 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, lb_exec_space, &
633 0 : int_exec_space, batch_size, xc_minus)
634 : END IF
635 :
636 0 : numerical_trace = (xc_plus - xc_minus)/(2.0_dp*dx)/3.0_dp
637 0 : diff_trace = analytic_trace - numerical_trace
638 :
639 0 : iw = cp_logger_get_default_io_unit()
640 0 : IF (iw > 0) THEN
641 : WRITE (UNIT=iw, FMT="(/,T2,A,1X,ES11.4)") &
642 0 : "GAUXC| Molecular XC virial finite-difference dx", dx
643 : WRITE (UNIT=iw, FMT="(T2,A,3(1X,ES19.11))") &
644 0 : "GAUXC| Molecular XC virial FD 1/3 Trace", &
645 0 : analytic_trace, numerical_trace, diff_trace
646 : END IF
647 :
648 0 : DEALLOCATE (particle_set_minus, particle_set_plus)
649 :
650 : END SUBROUTINE debug_gauxc_molecular_virial
651 :
652 : ! **************************************************************************************************
653 : !> \brief prints a force-based molecular XC virial diagnostic from GauXC gradients
654 : !> \param exc_grad ...
655 : !> \param particle_set ...
656 : !> \param para_env ...
657 : ! **************************************************************************************************
658 0 : SUBROUTINE print_gauxc_molecular_virial(exc_grad, particle_set, para_env)
659 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: exc_grad
660 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
661 : TYPE(mp_para_env_type), POINTER :: para_env
662 :
663 : CHARACTER(len=1), DIMENSION(3), PARAMETER :: label = ["x", "y", "z"]
664 :
665 : INTEGER :: i, iatom, iw, j
666 : REAL(KIND=dp), DIMENSION(3) :: center, displacement, grad, grad_sum
667 : REAL(KIND=dp), DIMENSION(3, 3) :: molecular_virial
668 :
669 0 : CPASSERT(ASSOCIATED(particle_set))
670 0 : CPASSERT(SIZE(exc_grad) == 3*SIZE(particle_set))
671 :
672 0 : IF (para_env%mepos /= 0) RETURN
673 :
674 0 : center = 0.0_dp
675 0 : DO iatom = 1, SIZE(particle_set)
676 0 : center = center + particle_set(iatom)%r
677 : END DO
678 0 : center = center/REAL(SIZE(particle_set), dp)
679 :
680 0 : grad_sum = 0.0_dp
681 0 : molecular_virial = 0.0_dp
682 0 : DO iatom = 1, SIZE(particle_set)
683 0 : grad = exc_grad(3*iatom - 2:3*iatom)
684 0 : displacement = particle_set(iatom)%r - center
685 0 : grad_sum = grad_sum + grad
686 0 : DO i = 1, 3
687 0 : DO j = 1, 3
688 0 : molecular_virial(i, j) = molecular_virial(i, j) + grad(i)*displacement(j)
689 : END DO
690 : END DO
691 : END DO
692 :
693 0 : iw = cp_logger_get_default_io_unit()
694 0 : IF (iw <= 0) RETURN
695 :
696 : WRITE (UNIT=iw, FMT="(/,T2,A)") &
697 0 : "GAUXC| Molecular XC gradient virial diagnostic [a.u.]"
698 0 : WRITE (UNIT=iw, FMT="(T2,A,T20,A,T40,A,T60,A)") "GAUXC|", "x", "y", "z"
699 0 : DO i = 1, 3
700 : WRITE (UNIT=iw, FMT="(T2,A,1X,A1,3(1X,ES19.11))") &
701 0 : "GAUXC|", label(i), molecular_virial(i, :)
702 : END DO
703 : WRITE (UNIT=iw, FMT="(T2,A,1X,ES19.11)") &
704 0 : "GAUXC| Molecular XC gradient virial 1/3 Trace", &
705 0 : (molecular_virial(1, 1) + molecular_virial(2, 2) + molecular_virial(3, 3))/3.0_dp
706 : WRITE (UNIT=iw, FMT="(T2,A,3(1X,ES19.11))") &
707 0 : "GAUXC| Molecular XC gradient sum", grad_sum
708 : WRITE (UNIT=iw, FMT="(T2,A)") &
709 0 : "GAUXC| Diagnostic only; this is not an analytical periodic stress tensor."
710 :
711 : END SUBROUTINE print_gauxc_molecular_virial
712 :
713 : ! **************************************************************************************************
714 : !> \brief Return information about the Skala functional
715 : !> \param functional section containing the SKALA subsection
716 : !> \param lsd if you are using lsd or lda
717 : !> \param reference the reference to the article where the functional is explained
718 : !> \param shortform the short definition of the functional
719 : !> \param needs the flags corresponding to the inputs needed by this
720 : !> functional are set to true (the flags not needed aren't touched)
721 : !> \param max_deriv the maximal derivative available
722 : ! **************************************************************************************************
723 68 : SUBROUTINE skala_info(functional, lsd, reference, shortform, needs, max_deriv)
724 : TYPE(section_vals_type), POINTER :: functional
725 : LOGICAL, INTENT(in) :: lsd
726 : CHARACTER(LEN=*), INTENT(OUT), OPTIONAL :: reference, shortform
727 : TYPE(xc_rho_cflags_type), INTENT(inout), OPTIONAL :: needs
728 : INTEGER, INTENT(out), OPTIONAL :: max_deriv
729 :
730 : CHARACTER(len=default_path_length) :: model_key, model_name
731 : CHARACTER(len=default_string_length) :: xc_fun_name
732 :
733 68 : CALL section_vals_val_get(functional, "FUNCTIONAL", c_val=xc_fun_name)
734 68 : CALL section_vals_val_get(functional, "MODEL", c_val=model_name)
735 68 : model_key = ADJUSTL(model_name)
736 68 : CALL uppercase(model_key)
737 :
738 68 : IF (PRESENT(reference)) THEN
739 8 : IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "") THEN
740 1 : reference = "Functional computed by GauXC (underlying: "//TRIM(xc_fun_name)//")"
741 : ELSE
742 7 : reference = "Functional computed by GauXC OneDFT model "//TRIM(model_name)
743 : END IF
744 : END IF
745 68 : IF (PRESENT(shortform)) THEN
746 8 : IF (TRIM(model_key) == "NONE" .OR. TRIM(model_key) == "") THEN
747 1 : shortform = "GAUXC ("//TRIM(xc_fun_name)//")"
748 : ELSE
749 7 : shortform = "GAUXC OneDFT"
750 : END IF
751 : END IF
752 68 : IF (PRESENT(needs)) THEN
753 60 : needs%rho = .TRUE.
754 60 : IF (lsd) THEN
755 20 : needs%rho_spin = .TRUE.
756 : END IF
757 : END IF
758 68 : IF (PRESENT(max_deriv)) max_deriv = 1
759 :
760 68 : END SUBROUTINE skala_info
761 :
762 : ! GauXC uses replicated dense density and VXC matrices. The DBCSR density matrix
763 : ! is distributed over MPI ranks, so apply_gauxc allreduces the dense copy before
764 : ! passing it to GauXC.
765 :
766 : ! **************************************************************************************************
767 : !> \brief ...
768 : !> \param qs_env ...
769 : !> \param xc_section ...
770 : !> \param calculate_forces ...
771 : ! **************************************************************************************************
772 510 : SUBROUTINE apply_gauxc(qs_env, xc_section, calculate_forces)
773 : TYPE(qs_environment_type), INTENT(in), POINTER :: qs_env
774 : TYPE(section_vals_type), INTENT(in), POINTER :: xc_section
775 : LOGICAL, INTENT(IN) :: calculate_forces
776 :
777 : CHARACTER(len=*), PARAMETER :: gapw_xc_abort_message = &
778 : "GauXC with METHOD GAPW_XC is not supported yet. "// &
779 : "The GAPW_XC one-center XC correction needs a dedicated GauXC design.", &
780 : nonlocal_vdw_abort_message = &
781 : "GauXC does not support non-local VDW_POTENTIAL corrections. "// &
782 : "Use an additive PAIR_POTENTIAL dispersion correction or disable GauXC."
783 : REAL(KIND=dp), PARAMETER :: gapw_fd_gradient_dx = 1.0E-4_dp
784 :
785 : CHARACTER(len=default_path_length) :: model_key, model_name, output_path
786 : CHARACTER(len=default_string_length) :: grid_key, grid_type, int_exec_space, lb_exec_space, &
787 : pruning_key, pruning_scheme, radial_quadrature, xc_fun_name
788 : INTEGER :: batch_size, env_status, img, ispin, &
789 : natom, nimages, nspins
790 : LOGICAL :: grid_explicit, hdf5_output, molecular_virial, molecular_virial_debug, &
791 : pruning_explicit, use_fd_gradient, use_gradient_self_runtime, use_onedft, &
792 : use_self_runtime, use_skala_model, write_hdf5_output
793 : REAL(KIND=dp) :: molecular_virial_debug_dx
794 510 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: density_scalar, density_zeta
795 510 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
796 : TYPE(cp_gauxc_basisset_type) :: gauxc_basis
797 : TYPE(cp_gauxc_grid_type) :: gauxc_gradient_grid_result, &
798 : gauxc_grid_result
799 : TYPE(cp_gauxc_integrator_type) :: gauxc_gradient_integrator_result, gauxc_integrator_result
800 : TYPE(cp_gauxc_molecule_type) :: gauxc_mol
801 : TYPE(cp_gauxc_status_type) :: gauxc_status
802 510 : TYPE(cp_gauxc_xc_gradient_type) :: exc_grad
803 510 : TYPE(cp_gauxc_xc_type) :: gauxc_xc_result
804 : TYPE(dbcsr_p_type) :: vxc_zeta_tmp
805 510 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vxc
806 510 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
807 : TYPE(dft_control_type), POINTER :: dft_control
808 : TYPE(mp_para_env_type), POINTER :: para_env
809 510 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
810 : TYPE(qs_energy_type), POINTER :: energy
811 510 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
812 510 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
813 : TYPE(qs_ks_env_type), POINTER :: ks_env
814 : TYPE(qs_rho_type), POINTER :: rho, rho_use, rho_xc
815 : TYPE(qs_scf_env_type), POINTER :: scf_env
816 : TYPE(section_vals_type), POINTER :: gauxc_functional_section
817 :
818 : NULLIFY ( &
819 : atomic_kind_set, &
820 510 : dft_control, &
821 510 : energy, &
822 510 : force, &
823 510 : ks_env, &
824 510 : matrix_vxc, &
825 510 : para_env, &
826 510 : particle_set, &
827 510 : qs_kind_set, &
828 510 : rho, &
829 510 : rho_use, &
830 510 : rho_xc, &
831 510 : rho_ao, &
832 510 : scf_env)
833 :
834 : CALL get_qs_env( &
835 : qs_env, &
836 : dft_control=dft_control, &
837 : energy=energy, &
838 : ks_env=ks_env, &
839 : matrix_vxc=matrix_vxc, &
840 : natom=natom, &
841 : atomic_kind_set=atomic_kind_set, &
842 : force=force, &
843 : para_env=para_env, &
844 : particle_set=particle_set, &
845 : qs_kind_set=qs_kind_set, &
846 : rho=rho, &
847 : rho_xc=rho_xc, &
848 510 : scf_env=scf_env)
849 :
850 510 : IF (dft_control%qs_control%gapw_xc) THEN
851 0 : CPABORT(gapw_xc_abort_message)
852 : END IF
853 510 : IF (dft_control%qs_control%gapw) THEN
854 52 : CALL ensure_gauxc_gapw_all_electron(qs_kind_set)
855 : END IF
856 510 : CPASSERT(ASSOCIATED(rho))
857 510 : rho_use => rho
858 : CALL qs_rho_get( &
859 : rho_use, &
860 510 : rho_ao_kp=rho_ao)
861 :
862 510 : nimages = dft_control%nimages
863 510 : nspins = dft_control%nspins
864 :
865 510 : IF (ASSOCIATED(qs_env%dispersion_env)) THEN
866 510 : IF (qs_env%dispersion_env%type == xc_vdw_fun_nonloc) THEN
867 0 : CPABORT(nonlocal_vdw_abort_message)
868 : END IF
869 : END IF
870 510 : NULLIFY (vxc_zeta_tmp%matrix)
871 :
872 510 : gauxc_functional_section => get_gauxc_functional(xc_section)
873 : CALL section_vals_val_get( &
874 : gauxc_functional_section, &
875 : "FUNCTIONAL", &
876 510 : c_val=xc_fun_name)
877 : CALL section_vals_val_get( &
878 : gauxc_functional_section, &
879 : "MODEL", &
880 510 : c_val=model_name)
881 : CALL section_vals_val_get( &
882 : gauxc_functional_section, &
883 : "GRID", &
884 : c_val=grid_type, &
885 510 : explicit=grid_explicit)
886 : CALL section_vals_val_get( &
887 : gauxc_functional_section, &
888 : "RADIAL_QUADRATURE", &
889 510 : c_val=radial_quadrature)
890 : CALL section_vals_val_get( &
891 : gauxc_functional_section, &
892 : "PRUNING_SCHEME", &
893 : c_val=pruning_scheme, &
894 510 : explicit=pruning_explicit)
895 : CALL section_vals_val_get( &
896 : gauxc_functional_section, &
897 : "BATCH_SIZE", &
898 510 : i_val=batch_size)
899 : CALL section_vals_val_get( &
900 : gauxc_functional_section, &
901 : "MOLECULAR_VIRIAL", &
902 510 : l_val=molecular_virial)
903 : CALL section_vals_val_get( &
904 : gauxc_functional_section, &
905 : "MOLECULAR_VIRIAL_DEBUG", &
906 510 : l_val=molecular_virial_debug)
907 : CALL section_vals_val_get( &
908 : gauxc_functional_section, &
909 : "MOLECULAR_VIRIAL_DEBUG_DX", &
910 510 : r_val=molecular_virial_debug_dx)
911 : CALL section_vals_val_get( &
912 : gauxc_functional_section, &
913 : "LB_EXECUTION_SPACE", &
914 510 : c_val=lb_exec_space)
915 : CALL section_vals_val_get( &
916 : gauxc_functional_section, &
917 : "INT_EXECUTION_SPACE", &
918 510 : c_val=int_exec_space)
919 : CALL section_vals_val_get( &
920 : gauxc_functional_section, &
921 : "OUTPUT_PATH", &
922 510 : c_val=output_path)
923 :
924 510 : model_key = ADJUSTL(model_name)
925 510 : CALL uppercase(model_key)
926 510 : use_onedft = (TRIM(model_key) /= "" .AND. TRIM(model_key) /= "NONE")
927 510 : use_skala_model = (INDEX(TRIM(model_key), "SKALA") > 0)
928 510 : IF (molecular_virial_debug) THEN
929 0 : IF (molecular_virial_debug_dx <= 0.0_dp) THEN
930 : CALL cp_abort(__LOCATION__, &
931 0 : "GauXC MOLECULAR_VIRIAL_DEBUG_DX must be positive.")
932 : END IF
933 0 : molecular_virial = .TRUE.
934 : END IF
935 510 : IF (dft_control%qs_control%gapw_xc .AND. use_onedft) THEN
936 : CALL cp_abort(__LOCATION__, &
937 : "GauXC OneDFT/SKALA with METHOD GAPW_XC is not implemented. "// &
938 0 : "The GAPW_XC one-center XC correction must be evaluated by the same GauXC model.")
939 : END IF
940 510 : IF (use_onedft) THEN
941 486 : IF (.NOT. grid_explicit) grid_type = "SUPERFINE"
942 486 : IF (.NOT. pruning_explicit) pruning_scheme = "UNPRUNED"
943 :
944 486 : grid_key = ADJUSTL(grid_type)
945 486 : pruning_key = ADJUSTL(pruning_scheme)
946 486 : CALL uppercase(grid_key)
947 486 : CALL uppercase(pruning_key)
948 486 : IF (use_skala_model .AND. (calculate_forces .OR. molecular_virial) .AND. &
949 : (TRIM(grid_key) /= "SUPERFINE" .OR. TRIM(pruning_key) /= "UNPRUNED")) THEN
950 : CALL cp_warn( &
951 : __LOCATION__, &
952 : "GauXC OneDFT/SKALA nuclear gradients are sensitive to the GauXC molecular grid. "// &
953 0 : "Use GRID SUPERFINE and PRUNING_SCHEME UNPRUNED for quantitative force checks.")
954 : END IF
955 486 : IF (TRIM(model_key) == "SKALA") THEN
956 152 : CALL GET_ENVIRONMENT_VARIABLE("GAUXC_SKALA_MODEL", model_name, STATUS=env_status)
957 152 : IF (env_status /= 0 .OR. LEN_TRIM(model_name) == 0) THEN
958 0 : CPABORT("MODEL SKALA requires the GAUXC_SKALA_MODEL environment variable")
959 : END IF
960 : END IF
961 : END IF
962 510 : use_self_runtime = use_skala_model .AND. para_env%num_pe > 1
963 : use_gradient_self_runtime = calculate_forces .AND. use_onedft .AND. &
964 510 : para_env%num_pe > 1 .AND. .NOT. use_self_runtime
965 :
966 : gauxc_mol = gauxc_create_molecule( &
967 : particle_set, &
968 510 : gauxc_status)
969 510 : CALL gauxc_check_status(gauxc_status)
970 : gauxc_basis = gauxc_create_basisset( &
971 : qs_kind_set, &
972 : particle_set, &
973 510 : gauxc_status)
974 510 : CALL gauxc_check_status(gauxc_status)
975 510 : hdf5_output = (TRIM(output_path) /= "")
976 510 : write_hdf5_output = hdf5_output .AND. para_env%mepos == 0
977 0 : IF (write_hdf5_output .AND. ASSOCIATED(scf_env)) THEN
978 0 : write_hdf5_output = scf_env%iter_count == 1
979 : END IF
980 0 : IF (write_hdf5_output) THEN
981 : CALL gauxc_write_molecule_hdf5( &
982 : gauxc_mol, &
983 : output_path, &
984 : "molecule.h5", &
985 : "molecule", &
986 0 : gauxc_status)
987 0 : CALL gauxc_check_status(gauxc_status)
988 : CALL gauxc_write_basisset_hdf5( &
989 : gauxc_basis, &
990 : output_path, &
991 : "basisset.h5", &
992 : "basisset", &
993 0 : gauxc_status)
994 0 : CALL gauxc_check_status(gauxc_status)
995 : END IF
996 : use_fd_gradient = dft_control%qs_control%gapw .AND. calculate_forces .AND. &
997 510 : (use_skala_model .OR. gauxc_basis%max_l > 3)
998 0 : IF (use_fd_gradient .AND. para_env%mepos == 0) THEN
999 : CALL cp_warn( &
1000 : __LOCATION__, &
1001 : "Using finite-difference GauXC XC gradients for METHOD GAPW with SKALA or "// &
1002 : "all-electron basis functions beyond f shells. The upstream analytical GauXC "// &
1003 0 : "gradient path is not yet reliable for this case.")
1004 : END IF
1005 510 : IF (use_self_runtime) THEN
1006 : ! SKALA currently needs a replicated molecular runtime for reproducible
1007 : ! open-shell densities across CP2K MPI ranks.
1008 : gauxc_grid_result = gauxc_create_grid( &
1009 : gauxc_mol, &
1010 : gauxc_basis, &
1011 : grid_type, &
1012 : radial_quadrature, &
1013 : pruning_scheme, &
1014 : lb_exec_space, &
1015 : batch_size, &
1016 : gauxc_status, &
1017 152 : mpi_comm=mp_comm_self%get_handle())
1018 : ELSE
1019 : ! Use the QS force-evaluation communicator rather than GauXC's global
1020 : ! runtime. In mixed CDFT the diabatic states can be built on disjoint
1021 : ! MPI subgroups; using the global communicator would make GauXC wait
1022 : ! for ranks that are working on another state.
1023 : gauxc_grid_result = gauxc_create_grid( &
1024 : gauxc_mol, &
1025 : gauxc_basis, &
1026 : grid_type, &
1027 : radial_quadrature, &
1028 : pruning_scheme, &
1029 : lb_exec_space, &
1030 : batch_size, &
1031 : gauxc_status, &
1032 358 : mpi_comm=para_env%get_handle())
1033 : END IF
1034 510 : CALL gauxc_check_status(gauxc_status)
1035 : gauxc_integrator_result = gauxc_create_integrator( &
1036 : TRIM(xc_fun_name), &
1037 : gauxc_grid_result, &
1038 : int_exec_space, &
1039 : nspins, &
1040 510 : gauxc_status)
1041 510 : CALL gauxc_check_status(gauxc_status)
1042 :
1043 510 : IF (use_gradient_self_runtime) THEN
1044 : ! Upstream GauXC does not yet support OneDFT/SKALA nuclear gradients
1045 : ! on an MPI runtime. Keep the energy/VXC path on the normal MPI
1046 : ! runtime and use an isolated runtime only for the replicated gradient.
1047 : gauxc_gradient_grid_result = gauxc_create_grid( &
1048 : gauxc_mol, &
1049 : gauxc_basis, &
1050 : grid_type, &
1051 : radial_quadrature, &
1052 : pruning_scheme, &
1053 : lb_exec_space, &
1054 : batch_size, &
1055 : gauxc_status, &
1056 4 : mpi_comm=mp_comm_self%get_handle())
1057 4 : CALL gauxc_check_status(gauxc_status)
1058 : gauxc_gradient_integrator_result = gauxc_create_integrator( &
1059 : TRIM(xc_fun_name), &
1060 : gauxc_gradient_grid_result, &
1061 : int_exec_space, &
1062 : nspins, &
1063 4 : gauxc_status)
1064 4 : CALL gauxc_check_status(gauxc_status)
1065 : END IF
1066 :
1067 510 : IF (qs_env%run_rtp) THEN
1068 0 : CPABORT("GAUXC XC energy currently does not support real-time propagation")
1069 : END IF
1070 :
1071 510 : energy%exc = 0
1072 :
1073 510 : IF (ASSOCIATED(matrix_vxc)) CALL dbcsr_deallocate_matrix_set(matrix_vxc)
1074 510 : CALL dbcsr_allocate_matrix_set(matrix_vxc, nspins)
1075 :
1076 1020 : DO img = 1, nimages
1077 510 : IF (img > 1) THEN
1078 0 : CPABORT("UNIMPLEMENTED: Handling nimg>1 in k-point integration")
1079 : END IF
1080 510 : CALL dbcsr_to_dense(rho_ao(1, img), density_scalar)
1081 510 : CALL para_env%sum(density_scalar)
1082 510 : IF (nspins == 1) THEN
1083 : gauxc_xc_result = gauxc_compute_xc( &
1084 : gauxc_integrator_result, &
1085 : density_scalar, &
1086 : nspins=nspins, &
1087 : status=gauxc_status, &
1088 384 : model=TRIM(model_name))
1089 384 : CALL gauxc_check_status(gauxc_status)
1090 384 : IF (calculate_forces) THEN
1091 4 : IF (use_fd_gradient) THEN
1092 : CALL gauxc_xc_gradient_fd( &
1093 : particle_set, qs_kind_set, density_scalar, nspins, model_name, &
1094 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
1095 : lb_exec_space, int_exec_space, batch_size, gapw_fd_gradient_dx, &
1096 0 : para_env, exc_grad%exc_grad)
1097 4 : ELSE IF (use_gradient_self_runtime) THEN
1098 : exc_grad = gauxc_compute_xc_gradient( &
1099 : gauxc_gradient_integrator_result, &
1100 : density_scalar, &
1101 : nspins=nspins, &
1102 : natom=natom, &
1103 : status=gauxc_status, &
1104 4 : model=TRIM(model_name))
1105 : ELSE
1106 : exc_grad = gauxc_compute_xc_gradient( &
1107 : gauxc_integrator_result, &
1108 : density_scalar, &
1109 : nspins=nspins, &
1110 : natom=natom, &
1111 : status=gauxc_status, &
1112 0 : model=TRIM(model_name))
1113 : END IF
1114 4 : CALL gauxc_check_status(gauxc_status)
1115 : CALL add_gauxc_gradient_to_force( &
1116 : exc_grad%exc_grad, &
1117 : force, &
1118 : atomic_kind_set, &
1119 4 : para_env)
1120 4 : IF (molecular_virial) THEN
1121 0 : CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1122 : END IF
1123 4 : IF (molecular_virial_debug) THEN
1124 : CALL debug_gauxc_molecular_virial( &
1125 : exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, nspins, &
1126 : model_name, xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
1127 0 : lb_exec_space, int_exec_space, batch_size, molecular_virial_debug_dx, para_env)
1128 : END IF
1129 4 : DEALLOCATE (exc_grad%exc_grad)
1130 : END IF
1131 : ELSE
1132 126 : CPASSERT(nspins == 2)
1133 : ! In here:
1134 : ! scalar <- rho_ao(1, :) + rho_ao(2, :)
1135 : ! zeta <- rho_ao(1, :) - rho_ao(2, :)
1136 126 : CALL dbcsr_to_dense(rho_ao(2, img), density_zeta)
1137 126 : CALL para_env%sum(density_zeta)
1138 : ! Do NOT reorder the following lines!
1139 26666 : density_scalar(:, :) = density_scalar(:, :) + density_zeta(:, :)
1140 : ! Factor two because the next line is evaluated after the above line.
1141 : ! We need to subtract density_zeta once to undo the above line and
1142 : ! a second time because that is what UKS requires.
1143 : ! This style lowers memory footprint.
1144 26666 : density_zeta(:, :) = density_scalar(:, :) - 2.0_dp*density_zeta(:, :)
1145 : gauxc_xc_result = gauxc_compute_xc( &
1146 : gauxc_integrator_result, &
1147 : density_scalar, &
1148 : density_zeta, &
1149 : nspins, &
1150 : gauxc_status, &
1151 126 : model=TRIM(model_name))
1152 126 : CALL gauxc_check_status(gauxc_status)
1153 126 : IF (calculate_forces) THEN
1154 0 : IF (use_fd_gradient) THEN
1155 : CALL gauxc_xc_gradient_fd( &
1156 : particle_set, qs_kind_set, density_scalar, nspins, model_name, &
1157 : xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
1158 : lb_exec_space, int_exec_space, batch_size, gapw_fd_gradient_dx, &
1159 0 : para_env, exc_grad%exc_grad, density_zeta=density_zeta)
1160 0 : ELSE IF (use_gradient_self_runtime) THEN
1161 : exc_grad = gauxc_compute_xc_gradient( &
1162 : gauxc_gradient_integrator_result, &
1163 : density_scalar, &
1164 : density_zeta, &
1165 : nspins, &
1166 : natom, &
1167 : gauxc_status, &
1168 0 : model=TRIM(model_name))
1169 : ELSE
1170 : exc_grad = gauxc_compute_xc_gradient( &
1171 : gauxc_integrator_result, &
1172 : density_scalar, &
1173 : density_zeta, &
1174 : nspins, &
1175 : natom, &
1176 : gauxc_status, &
1177 0 : model=TRIM(model_name))
1178 : END IF
1179 0 : CALL gauxc_check_status(gauxc_status)
1180 : CALL add_gauxc_gradient_to_force( &
1181 : exc_grad%exc_grad, &
1182 : force, &
1183 : atomic_kind_set, &
1184 0 : para_env)
1185 0 : IF (molecular_virial) THEN
1186 0 : CALL print_gauxc_molecular_virial(exc_grad%exc_grad, particle_set, para_env)
1187 : END IF
1188 0 : IF (molecular_virial_debug) THEN
1189 : CALL debug_gauxc_molecular_virial( &
1190 : exc_grad%exc_grad, particle_set, qs_kind_set, density_scalar, nspins, &
1191 : model_name, xc_fun_name, grid_type, radial_quadrature, pruning_scheme, &
1192 : lb_exec_space, int_exec_space, batch_size, molecular_virial_debug_dx, para_env, &
1193 0 : density_zeta=density_zeta)
1194 : END IF
1195 0 : DEALLOCATE (exc_grad%exc_grad)
1196 : END IF
1197 : END IF
1198 :
1199 510 : energy%exc = energy%exc + gauxc_xc_result%exc
1200 :
1201 1020 : IF (nspins == 1) THEN
1202 384 : IF (img == 1) THEN
1203 384 : matrix_vxc(1) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(1, img))
1204 : ELSE
1205 0 : CPABORT("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1206 : END IF
1207 : ELSE
1208 126 : CPASSERT(nspins == 2)
1209 : ! Transform derivatives from total/spin density back to alpha/beta channels.
1210 126 : vxc_zeta_tmp = dense_to_dbcsr(gauxc_xc_result%vxc_zeta, rho_ao(1, img))
1211 126 : IF (img == 1) THEN
1212 378 : DO ispin = 1, 2
1213 252 : matrix_vxc(ispin) = dense_to_dbcsr(gauxc_xc_result%vxc_scalar, rho_ao(ispin, 1))
1214 : CALL dbcsr_add( &
1215 : matrix_vxc(ispin)%matrix, &
1216 : vxc_zeta_tmp%matrix, &
1217 : 1.0_dp, &
1218 : ! 1.0 for ispin==1, -1.0 for ispin==2
1219 378 : 1.0_dp - REAL(ispin - 1, dp)*2.0_dp)
1220 : END DO
1221 : ELSE
1222 0 : CPABORT("UNIMPLEMENTED: Handling multiple result matrices in k-point integration")
1223 : END IF
1224 126 : CALL dbcsr_release(vxc_zeta_tmp%matrix)
1225 126 : DEALLOCATE (vxc_zeta_tmp%matrix)
1226 : END IF
1227 : END DO
1228 :
1229 510 : DEALLOCATE (density_scalar)
1230 510 : IF (ALLOCATED(density_zeta)) DEALLOCATE (density_zeta)
1231 510 : DEALLOCATE (gauxc_xc_result%vxc_scalar)
1232 510 : IF (ALLOCATED(gauxc_xc_result%vxc_zeta)) DEALLOCATE (gauxc_xc_result%vxc_zeta)
1233 :
1234 510 : CALL set_ks_env(ks_env, matrix_vxc=matrix_vxc)
1235 1146 : DO ispin = 1, nspins
1236 1146 : CALL dbcsr_finalize(matrix_vxc(ispin)%matrix)
1237 : END DO
1238 :
1239 510 : IF (use_gradient_self_runtime) THEN
1240 4 : CALL gauxc_destroy_integrator(gauxc_gradient_integrator_result, gauxc_status)
1241 4 : CALL gauxc_check_status(gauxc_status)
1242 4 : CALL gauxc_destroy_grid(gauxc_gradient_grid_result, gauxc_status)
1243 4 : CALL gauxc_check_status(gauxc_status)
1244 : END IF
1245 510 : CALL gauxc_destroy_integrator(gauxc_integrator_result, gauxc_status)
1246 510 : CALL gauxc_check_status(gauxc_status)
1247 510 : CALL gauxc_destroy_grid(gauxc_grid_result, gauxc_status)
1248 510 : CALL gauxc_check_status(gauxc_status)
1249 510 : CALL gauxc_destroy_basisset(gauxc_basis, gauxc_status)
1250 510 : CALL gauxc_check_status(gauxc_status)
1251 510 : CALL gauxc_destroy_molecule(gauxc_mol, gauxc_status)
1252 510 : CALL gauxc_check_status(gauxc_status)
1253 :
1254 1020 : END SUBROUTINE apply_gauxc
1255 :
1256 : END MODULE xc_gauxc_functional
|
