(git:6768d83)
Loading...
Searching...
No Matches
xc_gauxc_interface.F
Go to the documentation of this file.
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#define GAUXC_RETURN_IF_ERROR(status) IF (status%status%code /= 0) RETURN
13
14MODULE xc_gauxc_interface
15
16 USE iso_c_binding, ONLY: &
17 c_associated, &
18 c_bool, &
19 c_char, &
20 c_double, &
21 c_f_pointer, &
22 c_int, &
23 c_int32_t, &
24 c_int64_t, &
25 c_null_char, &
26 c_null_ptr, &
27 c_ptr, &
28 c_size_t
29 USE particle_types, ONLY: &
30 particle_type
31 USE qs_kind_types, ONLY: &
32 get_qs_kind, &
33 qs_kind_type
34 USE cp_dbcsr_api, ONLY: &
35 dbcsr_p_type
36
37#ifdef __GAUXC
38
39 USE kinds, ONLY: &
40 default_path_length, &
41 default_string_length, &
42 dp
43 USE physcon, ONLY: &
44 bohr
45 USE atomic_kind_types, ONLY: &
46 atomic_kind_type, &
47 get_atomic_kind, &
48 get_atomic_kind_set
49 USE qs_integral_utils, ONLY: &
50 basis_set_list_setup
51 USE basis_set_types, ONLY: &
52 gto_basis_set_p_type, &
53 gto_basis_set_type, &
54 write_gto_basis_set
55 USE periodic_table, ONLY: &
56 get_ptable_info
57 USE gauxc_status, ONLY: &
58 gauxc_status_message, &
59 gauxc_status_type
60 USE gauxc_enums, ONLY: &
61 gauxc_atomicgridsizedefault, &
62 gauxc_executionspace, &
63 gauxc_pruningscheme, &
64 gauxc_radialquad
65 USE gauxc_runtime_environment, ONLY: &
66 gauxc_runtime_environment_delete, &
67 gauxc_runtime_environment_new, &
68 gauxc_runtime_environment_type
69 USE gauxc_molecule, ONLY: &
70 gauxc_delete, &
71 gauxc_molecule_new_from_atoms, &
72 gauxc_molecule_type
73 USE gauxc_atom, ONLY: &
74 gauxc_atom_type
75 USE gauxc_basisset, ONLY: &
76 gauxc_basisset_new, &
77 gauxc_basisset_new_from_shells, &
78 gauxc_basisset_type, &
79 gauxc_delete
80 USE gauxc_shell, ONLY: &
81 gauxc_shell_type
82 USE gauxc_molgrid, ONLY: &
83 gauxc_delete, &
84 gauxc_molgrid_new_default, &
85 gauxc_molgrid_type
86 USE gauxc_load_balancer, ONLY: &
87 gauxc_delete, &
88 gauxc_load_balancer_factory_get_instance, &
89 gauxc_load_balancer_factory_new, &
90 gauxc_load_balancer_factory_type, &
91 gauxc_load_balancer_type
92 USE gauxc_molecular_weights, ONLY: &
93 gauxc_delete, &
94 gauxc_get_instance, &
95 gauxc_molecular_weights_factory_new, &
96 gauxc_molecular_weights_factory_type, &
97 gauxc_molecular_weights_modify_weights, &
98 gauxc_molecular_weights_settings, &
99 gauxc_molecular_weights_type
100 USE gauxc_xc_functional, ONLY: &
101 gauxc_delete, &
102 gauxc_functional_from_string, &
103 gauxc_functional_type
104 USE gauxc_integrator, ONLY: &
105 gauxc_delete, &
106 gauxc_integrator_eval_exc_grad_rks, &
107 gauxc_integrator_eval_exc_grad_uks, &
108 gauxc_integrator_eval_exc_vxc_rks, &
109 gauxc_integrator_eval_exc_vxc_uks, &
110 gauxc_integrator_new, &
111 gauxc_integrator_type
112#ifdef GAUXC_HAS_ONEDFT
113 USE gauxc_integrator, ONLY: &
114 gauxc_integrator_eval_exc_grad_onedft_uks, &
115 gauxc_integrator_eval_exc_vxc_onedft_uks
116 USE omp_lib, ONLY: &
117 omp_get_max_threads, &
118 omp_set_num_threads
119#endif
120 USE string_utilities, ONLY: &
121 uppercase
122#endif
123
124#include "../base/base_uses.f90"
125
126 IMPLICIT NONE
127 PRIVATE
128
129#ifndef __GAUXC
130
131 ! The module still exists as an empty shell when compiling without GauXC.
132
133 TYPE cp_gauxc_molecule_type
134 END TYPE cp_gauxc_molecule_type
135
136 TYPE cp_gauxc_basisset_type
137 INTEGER :: max_l = -1
138 END TYPE cp_gauxc_basisset_type
139
140 TYPE cp_gauxc_grid_type
141 END TYPE cp_gauxc_grid_type
142
143 TYPE cp_gauxc_integrator_type
144 END TYPE cp_gauxc_integrator_type
145
146 TYPE cp_gauxc_status_type
147 END TYPE cp_gauxc_status_type
148
149#else
150
151 ! TODO can we make the single fields private somehow?
152
153 TYPE cp_gauxc_molecule_type
154 TYPE(gauxc_molecule_type) :: molecule
155 END TYPE cp_gauxc_molecule_type
156
157 TYPE cp_gauxc_basisset_type
158 TYPE(gauxc_basisset_type) :: basis
159 INTEGER :: max_l = -1
160 END TYPE cp_gauxc_basisset_type
161
162 TYPE cp_gauxc_grid_type
163 TYPE(gauxc_molgrid_type) :: grid
164 TYPE(gauxc_load_balancer_type) :: lb
165 TYPE(gauxc_load_balancer_factory_type) :: lbf
166 TYPE(gauxc_molecular_weights_type) :: mw
167 TYPE(gauxc_molecular_weights_factory_type) :: mwf
168 TYPE(gauxc_runtime_environment_type) :: rt
169 LOGICAL :: owns_rt = .false.
170 END TYPE cp_gauxc_grid_type
171
172 TYPE cp_gauxc_integrator_type
173 TYPE(gauxc_functional_type) :: func
174 TYPE(gauxc_integrator_type) :: integrator
175 END TYPE cp_gauxc_integrator_type
176
177 TYPE cp_gauxc_status_type
178 TYPE(gauxc_status_type) :: status
179 END TYPE cp_gauxc_status_type
180
181 TYPE(gauxc_runtime_environment_type) :: rt
182 INTEGER :: rt_mpi_comm = -1
183 LOGICAL :: rt_has_mpi_comm = .false.
184
185#endif
186
187 TYPE cp_gauxc_xc_type
188 REAL(c_double) :: exc
189 REAL(c_double), DIMENSION(:, :), ALLOCATABLE :: vxc_scalar, vxc_zeta
190 END TYPE cp_gauxc_xc_type
191
192 TYPE cp_gauxc_xc_gradient_type
193 REAL(c_double), ALLOCATABLE, DIMENSION(:) :: exc_grad
194 END TYPE cp_gauxc_xc_gradient_type
195
196 CHARACTER(len=*), PARAMETER :: no_gauxc_message = "Compile CP2K with GauXC to use this functionality!"
197
198 PUBLIC :: &
199 cp_gauxc_basisset_type, &
200 cp_gauxc_grid_type, &
201 cp_gauxc_integrator_type, &
202 cp_gauxc_molecule_type, &
203 cp_gauxc_status_type, &
204 cp_gauxc_xc_gradient_type, &
205 cp_gauxc_xc_type, &
206 gauxc_check_status, &
207 gauxc_compute_xc_gradient, &
208 gauxc_compute_xc, &
209 gauxc_create_basisset, &
210 gauxc_create_grid, &
211 gauxc_create_integrator, &
212 gauxc_create_molecule, &
213 gauxc_destroy_basisset, &
214 gauxc_destroy_grid, &
215 gauxc_destroy_integrator, &
216 gauxc_destroy_molecule, &
217 gauxc_finalize, &
218 gauxc_init
219CONTAINS
220
221! **************************************************************************************************
222!> \brief ...
223!> \param status ...
224! **************************************************************************************************
225 SUBROUTINE print_gauxc_status_message(status)
226 ! IMPORT :: c_ptr
227 TYPE(cp_gauxc_status_type) :: status
228
229#ifdef __GAUXC
230 CHARACTER(kind=c_char), POINTER :: s(:)
231 INTEGER :: i
232
233 WRITE (0, '(a,1x,i0)') "GauXC returned with status code", status%status%code
234 IF (c_associated(status%status%message)) THEN
235 WRITE (0, '(a)', advance='no') "GauXC status message: ["
236
237 CALL c_f_pointer(status%status%message, s, [default_string_length])
238 DO i = 1, SIZE(s)
239 IF (s(i) == c_null_char) EXIT
240 WRITE (0, '(A)', advance='no') s(i)
241 END DO
242
243 WRITE (0, '(a)') "]"
244 ELSE
245 WRITE (0, '(a)') "GauXC status message: [null]"
246 END IF
247#else
248 mark_used(status)
249#endif
250 END SUBROUTINE print_gauxc_status_message
251
252! **************************************************************************************************
253!> \brief ...
254!> \param mpi_comm ...
255!> \param status ...
256! **************************************************************************************************
257 SUBROUTINE gauxc_init(mpi_comm, status)
258 INTEGER, INTENT(IN), OPTIONAL :: mpi_comm
259 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
260
261#ifdef __GAUXC
262#if defined(GAUXC_HAS_MPI) && defined(__parallel)
263 IF (PRESENT(mpi_comm)) THEN
264 rt = gauxc_runtime_environment_new(status%status, mpi_comm)
265 rt_mpi_comm = mpi_comm
266 rt_has_mpi_comm = .true.
267 ELSE
268 rt = gauxc_runtime_environment_new(status%status)
269 rt_mpi_comm = -1
270 rt_has_mpi_comm = .false.
271 END IF
272#else
273 mark_used(mpi_comm)
274 rt = gauxc_runtime_environment_new(status%status)
275 rt_mpi_comm = -1
276 rt_has_mpi_comm = .false.
277#endif
278 gauxc_return_if_error(status)
279#else
280 mark_used(mpi_comm)
281 mark_used(status)
282#endif
283 END SUBROUTINE gauxc_init
284
285! **************************************************************************************************
286!> \brief ...
287!> \param status ...
288! **************************************************************************************************
289 SUBROUTINE gauxc_finalize(status)
290 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
291
292#ifdef __GAUXC
293 CALL gauxc_runtime_environment_delete(status%status, rt)
294 gauxc_return_if_error(status)
295 rt_mpi_comm = -1
296 rt_has_mpi_comm = .false.
297#else
298 mark_used(status)
299#endif
300 END SUBROUTINE gauxc_finalize
301
302! **************************************************************************************************
303!> \brief ...
304!> \param particle_set ...
305!> \param status ...
306!> \return ...
307! **************************************************************************************************
308 FUNCTION gauxc_create_molecule(particle_set, status) RESULT(res)
309 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
310 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
311 TYPE(cp_gauxc_molecule_type) :: res
312
313#ifdef __GAUXC
314 CHARACTER(LEN=2) :: element_symbol
315 INTEGER :: atomic_number, i, natoms
316 TYPE(atomic_kind_type), POINTER :: atomic_kind
317 TYPE(gauxc_atom_type), ALLOCATABLE, DIMENSION(:) :: atoms
318
319 natoms = SIZE(particle_set)
320 ALLOCATE (atoms(natoms))
321
322 DO i = 1, natoms
323 atomic_kind => particle_set(i)%atomic_kind
324 CALL get_atomic_kind(atomic_kind, element_symbol=element_symbol)
325 CALL get_ptable_info(element_symbol, number=atomic_number)
326 atoms(i)%atomic_number = int(atomic_number, c_int64_t)
327 atoms(i)%x = real(particle_set(i)%r(1), c_double)
328 atoms(i)%y = real(particle_set(i)%r(2), c_double)
329 atoms(i)%z = real(particle_set(i)%r(3), c_double)
330 END DO
331
332 res%molecule = gauxc_molecule_new_from_atoms(status%status, atoms, int(natoms, c_size_t))
333 gauxc_return_if_error(status)
334
335 DEALLOCATE (atoms)
336#else
337 mark_used(particle_set)
338 mark_used(res)
339 mark_used(status)
340 cpabort(no_gauxc_message)
341#endif
342 END FUNCTION gauxc_create_molecule
343
344! **************************************************************************************************
345!> \brief ...
346!> \param qs_kind_set ...
347!> \param particle_set ...
348!> \param status ...
349!> \return ...
350! **************************************************************************************************
351 FUNCTION gauxc_create_basisset(qs_kind_set, particle_set, status) RESULT(res)
352 TYPE(qs_kind_type), DIMENSION(:), INTENT(IN), &
353 POINTER :: qs_kind_set
354 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
355 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
356 TYPE(cp_gauxc_basisset_type) :: res
357
358#ifdef __GAUXC
359 INTEGER :: iatom, ikind, iprim, iset, ishell, lval, &
360 nkind, npgf, nset, nshell, &
361 nshell_total, shell_index, natoms
362 REAL(c_double), DIMENSION(3) :: shell_origin
363 TYPE(atomic_kind_type), POINTER :: atomic_kind
364 TYPE(gauxc_shell_type), ALLOCATABLE, DIMENSION(:) :: shells
365 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
366 TYPE(gto_basis_set_type), POINTER :: gto_basis
367
368 nkind = SIZE(qs_kind_set)
369 natoms = SIZE(particle_set)
370
371 ALLOCATE (basis_set_list(nkind))
372 CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
373
374 nshell_total = 0
375 DO iatom = 1, natoms
376 atomic_kind => particle_set(iatom)%atomic_kind
377 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
378 gto_basis => basis_set_list(ikind)%gto_basis_set
379 cpassert(ASSOCIATED(gto_basis))
380 ! CALL write_gto_basis_set(gto_basis, 6, "–––––GAUXC-CREATE-BASISSET–––––")
381 nshell_total = nshell_total + sum(gto_basis%nshell)
382 END DO
383
384 ALLOCATE (shells(nshell_total))
385
386 shell_index = 0
387 res%max_l = -1
388 DO iatom = 1, natoms ! for each atom
389 atomic_kind => particle_set(iatom)%atomic_kind
390 CALL get_atomic_kind(atomic_kind, kind_number=ikind)
391 gto_basis => basis_set_list(ikind)%gto_basis_set
392 cpassert(ASSOCIATED(gto_basis))
393
394 shell_origin(1) = real(particle_set(iatom)%r(1), c_double)
395 shell_origin(2) = real(particle_set(iatom)%r(2), c_double)
396 shell_origin(3) = real(particle_set(iatom)%r(3), c_double)
397
398 nset = gto_basis%nset
399 cpassert(nset == SIZE(gto_basis%nshell))
400 DO iset = 1, nset ! for each shell group
401 nshell = gto_basis%nshell(iset)
402 npgf = gto_basis%npgf(iset) ! corresponds with nprim of gauxc
403
404 DO ishell = 1, gto_basis%nshell(iset) ! for each shell within the shell group
405 shell_index = shell_index + 1 ! global shell index, flattened over atoms and groups
406 lval = gto_basis%l(ishell, iset)
407 res%max_l = max(res%max_l, lval)
408 shells(shell_index)%l = int(lval, c_int32_t)
409 ! FIXME hardcoded true param
410 ! pure=1: spherical Gaussians; pure=0: cartesian Gaussians
411 shells(shell_index)%pure = .true._c_bool
412 shells(shell_index)%nprim = int(npgf, c_int32_t)
413 shells(shell_index)%origin = shell_origin
414
415 DO iprim = 1, npgf
416 shells(shell_index)%exponents(iprim) = &
417 REAL(gto_basis%zet(iprim, iset), c_double)
418 shells(shell_index)%coefficients(iprim) = &
419 REAL(gto_basis%norm_cgf(gto_basis%first_cgf(ishell, iset))* &
420 gto_basis%gcc(iprim, ishell, iset), c_double)
421 END DO
422 END DO
423 END DO
424 END DO
425
426 res%basis = gauxc_basisset_new_from_shells( &
427 status%status, &
428 shells, &
429 normalize=.false.)
430 gauxc_return_if_error(status)
431
432 DEALLOCATE (shells)
433 DEALLOCATE (basis_set_list)
434
435#else
436 mark_used(particle_set)
437 mark_used(qs_kind_set)
438 mark_used(res)
439 mark_used(status)
440 cpabort(no_gauxc_message)
441#endif
442 END FUNCTION gauxc_create_basisset
443
444! **************************************************************************************************
445!> \brief ...
446!> \param molecule ...
447!> \param basis ...
448!> \param grid_type ...
449!> \param radial_quadrature ...
450!> \param pruning_scheme ...
451!> \param lb_exec_space ...
452!> \param batch_size ...
453!> \param status ...
454!> \param mpi_comm optional communicator for a grid-local GauXC runtime
455!> \return ...
456! **************************************************************************************************
457 FUNCTION gauxc_create_grid( &
458 molecule, &
459 basis, &
460 grid_type, &
461 radial_quadrature, &
462 pruning_scheme, &
463 lb_exec_space, &
464 batch_size, &
465 status, &
466 mpi_comm) RESULT(res)
467
468 TYPE(cp_gauxc_molecule_type), INTENT(IN) :: molecule
469 TYPE(cp_gauxc_basisset_type), INTENT(in) :: basis
470 CHARACTER(len=*) :: grid_type, lb_exec_space, &
471 pruning_scheme, radial_quadrature
472 INTEGER :: batch_size
473 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
474 INTEGER, INTENT(IN), OPTIONAL :: mpi_comm
475 TYPE(cp_gauxc_grid_type) :: res
476
477#ifdef __GAUXC
478 INTEGER(c_int) :: grid_type_local, int_exec_space_local, &
479 lb_exec_space_local, &
480 pruning_scheme_local, radial_quad_local
481
482 grid_type_local = read_atomic_grid_size(grid_type)
483 radial_quad_local = read_radial_quad(radial_quadrature)
484 pruning_scheme_local = read_pruning_scheme(pruning_scheme)
485 lb_exec_space_local = read_execution_space(lb_exec_space)
486 int_exec_space_local = read_execution_space("host")
487 res%owns_rt = .false.
488
489#if defined(GAUXC_HAS_MPI) && defined(__parallel)
490 IF (PRESENT(mpi_comm)) THEN
491 ! Reuse the global runtime when the requested communicator matches
492 ! the communicator used during gauxc_init.
493 IF (.NOT. rt_has_mpi_comm .OR. mpi_comm /= rt_mpi_comm) THEN
494 res%rt = gauxc_runtime_environment_new(status%status, mpi_comm)
495 gauxc_return_if_error(status)
496 res%owns_rt = .true.
497 END IF
498 END IF
499#else
500 mark_used(mpi_comm)
501#endif
502
503 res%grid = gauxc_molgrid_new_default( &
504 status%status, &
505 molecule%molecule, &
506 pruning_scheme_local, &
507 int(batch_size, c_int64_t), &
508 radial_quad_local, &
509 grid_type_local)
510 gauxc_return_if_error(status)
511
512 res%lbf = gauxc_load_balancer_factory_new( &
513 status%status, &
514 lb_exec_space_local)
515 gauxc_return_if_error(status)
516
517 IF (res%owns_rt) THEN
518 res%lb = gauxc_load_balancer_factory_get_instance( &
519 status%status, &
520 res%lbf, &
521 res%rt, &
522 molecule%molecule, &
523 res%grid, &
524 basis%basis)
525 ELSE
526 res%lb = gauxc_load_balancer_factory_get_instance( &
527 status%status, &
528 res%lbf, &
529 rt, &
530 molecule%molecule, &
531 res%grid, &
532 basis%basis)
533 END IF
534 gauxc_return_if_error(status)
535
536 res%mwf = gauxc_molecular_weights_factory_new( &
537 status%status, &
538 int_exec_space_local)
539 gauxc_return_if_error(status)
540
541 res%mw = gauxc_get_instance( &
542 status%status, &
543 res%mwf)
544 gauxc_return_if_error(status)
545
546 CALL gauxc_molecular_weights_modify_weights( &
547 status%status, &
548 res%mw, &
549 res%lb)
550 gauxc_return_if_error(status)
551
552#else
553 mark_used(basis)
554 mark_used(batch_size)
555 mark_used(grid_type)
556 mark_used(lb_exec_space)
557 mark_used(mpi_comm)
558 mark_used(molecule)
559 mark_used(pruning_scheme)
560 mark_used(radial_quadrature)
561 mark_used(res)
562 mark_used(status)
563 cpabort(no_gauxc_message)
564#endif
565 END FUNCTION gauxc_create_grid
566
567! **************************************************************************************************
568!> \brief ...
569!> \param xc_functional_name ...
570!> \param grid ...
571!> \param int_exec_space ...
572!> \param nspins ...
573!> \param status ...
574!> \return ...
575! **************************************************************************************************
576 FUNCTION gauxc_create_integrator( &
577 xc_functional_name, &
578 grid, &
579 int_exec_space, &
580 nspins, &
581 status) RESULT(res)
582
583 CHARACTER(len=*), INTENT(IN) :: xc_functional_name, int_exec_space
584 TYPE(cp_gauxc_grid_type), INTENT(IN) :: grid
585 INTEGER, INTENT(IN) :: nspins
586 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
587 TYPE(cp_gauxc_integrator_type) :: res
588
589#ifdef __GAUXC
590 INTEGER(c_int) :: int_exec_space_local
591 LOGICAL(c_bool) :: polarized
592
593 polarized = (nspins == 2)
594 res%func = gauxc_functional_from_string( &
595 status%status, &
596 xc_functional_name, &
597 polarized)
598 gauxc_return_if_error(status)
599
600 int_exec_space_local = read_execution_space(int_exec_space)
601 res%integrator = gauxc_integrator_new( &
602 status%status, &
603 res%func, &
604 grid%lb, &
605 int_exec_space_local)
606 gauxc_return_if_error(status)
607
608#else
609 mark_used(grid)
610 mark_used(int_exec_space)
611 mark_used(nspins)
612 mark_used(res)
613 mark_used(status)
614 mark_used(xc_functional_name)
615 cpabort(no_gauxc_message)
616#endif
617 END FUNCTION gauxc_create_integrator
618
619! **************************************************************************************************
620!> \brief ...
621!> \param integrator ...
622!> \param density_scalar ...
623!> \param density_zeta ...
624!> \param nspins ...
625!> \param status ...
626!> \param model ...
627!> \return ...
628! **************************************************************************************************
629 FUNCTION gauxc_compute_xc( &
630 integrator, &
631 density_scalar, &
632 density_zeta, &
633 nspins, &
634 status, &
635 model) RESULT(res)
636
637 TYPE(cp_gauxc_integrator_type), INTENT(IN) :: integrator
638 REAL(c_double), DIMENSION(:, :), INTENT(IN) :: density_scalar
639 REAL(c_double), DIMENSION(:, :), INTENT(IN), &
640 OPTIONAL :: density_zeta
641 INTEGER, INTENT(IN) :: nspins
642 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
643 CHARACTER(len=*), INTENT(IN), OPTIONAL :: model
644 TYPE(cp_gauxc_xc_type) :: res
645
646#ifdef __GAUXC
647 CHARACTER(len=default_path_length) :: model_key
648 LOGICAL :: use_onedft
649#ifdef GAUXC_HAS_ONEDFT
650 REAL(c_double), ALLOCATABLE, DIMENSION(:, :) :: density_zeta_zero
651 INTEGER :: omp_max_threads_restore
652#endif
653
654 use_onedft = .false.
655 IF (PRESENT(model)) THEN
656 model_key = adjustl(model)
657 CALL uppercase(model_key)
658 use_onedft = (trim(model_key) /= "" .AND. trim(model_key) /= "NONE")
659 END IF
660
661 IF (.NOT. ALLOCATED(res%vxc_scalar)) THEN
662 ALLOCATE (res%vxc_scalar, mold=density_scalar)
663 ELSE
664 cpassert(all(shape(res%vxc_scalar) == shape(density_scalar)))
665 END IF
666 res%vxc_scalar = 0._dp
667
668 IF (use_onedft) THEN
669#ifndef GAUXC_HAS_ONEDFT
670 cpabort("GauXC lacks OneDFT support")
671#else
672 ! OneDFT may change the OpenMP team size for later parallel regions.
673 ! Restore max threads only; omp_get_num_threads() is 1 here.
674 omp_max_threads_restore = omp_get_max_threads()
675 IF (.NOT. ALLOCATED(res%vxc_zeta)) THEN
676 ALLOCATE (res%vxc_zeta, mold=density_scalar)
677 ELSE
678 cpassert(all(shape(res%vxc_zeta) == shape(density_scalar)))
679 END IF
680 res%vxc_zeta = 0._dp
681
682 IF (nspins == 1) THEN
683 ALLOCATE (density_zeta_zero, mold=density_scalar)
684 density_zeta_zero = 0._dp
685 CALL gauxc_integrator_eval_exc_vxc_onedft_uks( &
686 status%status, &
687 integrator%integrator, &
688 density_scalar, &
689 density_zeta_zero, &
690 trim(model), &
691 res%exc, &
692 res%vxc_scalar, &
693 res%vxc_zeta)
694 DEALLOCATE (density_zeta_zero)
695 ELSE
696 cpassert(PRESENT(density_zeta))
697 CALL gauxc_integrator_eval_exc_vxc_onedft_uks( &
698 status%status, &
699 integrator%integrator, &
700 density_scalar, &
701 density_zeta, &
702 trim(model), &
703 res%exc, &
704 res%vxc_scalar, &
705 res%vxc_zeta)
706 END IF
707 CALL omp_set_num_threads(omp_max_threads_restore)
708 gauxc_return_if_error(status)
709 RETURN
710#endif
711 END IF
712
713 IF (nspins == 1) THEN
714 CALL gauxc_integrator_eval_exc_vxc_rks( &
715 status%status, &
716 integrator%integrator, &
717 density_scalar, &
718 res%exc, &
719 res%vxc_scalar)
720 ELSE
721 cpassert(PRESENT(density_zeta))
722
723 IF (.NOT. ALLOCATED(res%vxc_zeta)) THEN
724 ALLOCATE (res%vxc_zeta, mold=density_zeta)
725 ELSE
726 cpassert(all(shape(res%vxc_zeta) == shape(density_scalar)))
727 END IF
728 res%vxc_zeta = 0._dp
729
730 CALL gauxc_integrator_eval_exc_vxc_uks( &
731 status%status, &
732 integrator%integrator, &
733 density_scalar, &
734 density_zeta, &
735 res%exc, &
736 res%vxc_scalar, &
737 res%vxc_zeta)
738 END IF
739 gauxc_return_if_error(status)
740
741#else
742 mark_used(integrator)
743 mark_used(density_scalar)
744 mark_used(density_zeta)
745 mark_used(nspins)
746 mark_used(status)
747 mark_used(model)
748 cpabort(no_gauxc_message)
749#endif
750 END FUNCTION gauxc_compute_xc
751
752! **************************************************************************************************
753!> \brief ...
754!> \param integrator ...
755!> \param density_scalar ...
756!> \param density_zeta ...
757!> \param nspins ...
758!> \param natom ...
759!> \param status ...
760!> \param model ...
761!> \return ...
762! **************************************************************************************************
763 FUNCTION gauxc_compute_xc_gradient( &
764 integrator, &
765 density_scalar, &
766 density_zeta, &
767 nspins, &
768 natom, &
769 status, &
770 model) RESULT(res)
771
772 TYPE(cp_gauxc_integrator_type), INTENT(IN) :: integrator
773 REAL(c_double), DIMENSION(:, :), INTENT(IN) :: density_scalar
774 REAL(c_double), DIMENSION(:, :), INTENT(IN), &
775 OPTIONAL :: density_zeta
776 INTEGER, INTENT(IN) :: nspins, natom
777 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
778 CHARACTER(len=*), INTENT(IN), OPTIONAL :: model
779 TYPE(cp_gauxc_xc_gradient_type) :: res
780
781#ifdef __GAUXC
782 CHARACTER(len=default_path_length) :: model_key
783 LOGICAL :: use_onedft
784#ifdef GAUXC_HAS_ONEDFT
785 REAL(c_double), ALLOCATABLE, DIMENSION(:, :) :: density_zeta_zero
786 INTEGER :: omp_max_threads_restore
787#endif
788
789 ALLOCATE (res%exc_grad(3*natom))
790 res%exc_grad = 0._dp
791
792 use_onedft = .false.
793 IF (PRESENT(model)) THEN
794 model_key = adjustl(model)
795 CALL uppercase(model_key)
796 use_onedft = (trim(model_key) /= "" .AND. trim(model_key) /= "NONE")
797 END IF
798
799 IF (use_onedft) THEN
800#ifndef GAUXC_HAS_ONEDFT
801 cpabort("GauXC lacks OneDFT support")
802#else
803 ! OneDFT may change the OpenMP team size for later parallel regions.
804 ! Restore max threads only; omp_get_num_threads() is 1 here.
805 omp_max_threads_restore = omp_get_max_threads()
806 IF (nspins == 1) THEN
807 ALLOCATE (density_zeta_zero, mold=density_scalar)
808 density_zeta_zero = 0._dp
809 CALL gauxc_integrator_eval_exc_grad_onedft_uks( &
810 status%status, &
811 integrator%integrator, &
812 density_scalar, &
813 density_zeta_zero, &
814 trim(model), &
815 res%exc_grad)
816 DEALLOCATE (density_zeta_zero)
817 ELSE
818 cpassert(PRESENT(density_zeta))
819 CALL gauxc_integrator_eval_exc_grad_onedft_uks( &
820 status%status, &
821 integrator%integrator, &
822 density_scalar, &
823 density_zeta, &
824 trim(model), &
825 res%exc_grad)
826 END IF
827 CALL omp_set_num_threads(omp_max_threads_restore)
828 gauxc_return_if_error(status)
829 RETURN
830#endif
831 END IF
832
833 IF (nspins == 1) THEN
834 CALL gauxc_integrator_eval_exc_grad_rks( &
835 status%status, &
836 integrator%integrator, &
837 density_scalar, &
838 res%exc_grad)
839 ELSE
840 cpassert(PRESENT(density_zeta))
841 CALL gauxc_integrator_eval_exc_grad_uks( &
842 status%status, &
843 integrator%integrator, &
844 density_scalar, &
845 density_zeta, &
846 res%exc_grad)
847 END IF
848 gauxc_return_if_error(status)
849
850#else
851 mark_used(density_scalar)
852 mark_used(density_zeta)
853 mark_used(res)
854 mark_used(integrator)
855 mark_used(model)
856 mark_used(natom)
857 mark_used(nspins)
858 mark_used(status)
859 cpabort(no_gauxc_message)
860#endif
861 END FUNCTION gauxc_compute_xc_gradient
862
863! **************************************************************************************************
864!> \brief ...
865!> \param molecule ...
866!> \param status ...
867! **************************************************************************************************
868 SUBROUTINE gauxc_destroy_molecule(molecule, status)
869 TYPE(cp_gauxc_molecule_type), INTENT(INOUT) :: molecule
870 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
871
872#ifdef __GAUXC
873 CALL gauxc_delete(status%status, molecule%molecule)
874 gauxc_return_if_error(status)
875#else
876 mark_used(molecule)
877 mark_used(status)
878 cpabort(no_gauxc_message)
879#endif
880 END SUBROUTINE gauxc_destroy_molecule
881
882! **************************************************************************************************
883!> \brief ...
884!> \param basis ...
885!> \param status ...
886! **************************************************************************************************
887 SUBROUTINE gauxc_destroy_basisset(basis, status)
888 TYPE(cp_gauxc_basisset_type), INTENT(INOUT) :: basis
889 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
890
891#ifdef __GAUXC
892 CALL gauxc_delete(status%status, basis%basis)
893 gauxc_return_if_error(status)
894#else
895 mark_used(basis)
896 mark_used(status)
897 cpabort(no_gauxc_message)
898#endif
899 END SUBROUTINE gauxc_destroy_basisset
900
901! **************************************************************************************************
902!> \brief ...
903!> \param grid_result ...
904!> \param status ...
905! **************************************************************************************************
906 SUBROUTINE gauxc_destroy_grid(grid_result, status)
907 TYPE(cp_gauxc_grid_type), INTENT(INOUT) :: grid_result
908 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
909
910#ifdef __GAUXC
911 CALL gauxc_delete(status%status, grid_result%mw)
912 gauxc_return_if_error(status)
913 CALL gauxc_delete(status%status, grid_result%mwf)
914 gauxc_return_if_error(status)
915 CALL gauxc_delete(status%status, grid_result%lb)
916 gauxc_return_if_error(status)
917 CALL gauxc_delete(status%status, grid_result%lbf)
918 gauxc_return_if_error(status)
919 CALL gauxc_delete(status%status, grid_result%grid)
920 gauxc_return_if_error(status)
921 IF (grid_result%owns_rt) THEN
922 CALL gauxc_runtime_environment_delete(status%status, grid_result%rt)
923 gauxc_return_if_error(status)
924 grid_result%owns_rt = .false.
925 END IF
926#else
927 mark_used(grid_result)
928 mark_used(status)
929 cpabort(no_gauxc_message)
930#endif
931 END SUBROUTINE gauxc_destroy_grid
932
933! **************************************************************************************************
934!> \brief ...
935!> \param integrator_result ...
936!> \param status ...
937! **************************************************************************************************
938 SUBROUTINE gauxc_destroy_integrator(integrator_result, status)
939 TYPE(cp_gauxc_integrator_type), INTENT(INOUT) :: integrator_result
940 TYPE(cp_gauxc_status_type), INTENT(OUT) :: status
941
942#ifdef __GAUXC
943 CALL gauxc_delete(status%status, integrator_result%integrator)
944 gauxc_return_if_error(status)
945 CALL gauxc_delete(status%status, integrator_result%func)
946 gauxc_return_if_error(status)
947#else
948 mark_used(integrator_result)
949 mark_used(status)
950 cpabort(no_gauxc_message)
951#endif
952 END SUBROUTINE gauxc_destroy_integrator
953
954! **************************************************************************************************
955!> \brief Checks gauxc status and prints error message before aborting
956!> \param status the gauxc status to check
957! **************************************************************************************************
958 SUBROUTINE gauxc_check_status(status)
959 TYPE(cp_gauxc_status_type), INTENT(IN) :: status
960
961#ifdef __GAUXC
962 IF (status%status%code /= 0) THEN
963 CALL print_gauxc_status_message(status)
964 cpabort("GauXC returned with non-zero status code")
965 END IF
966#else
967 mark_used(status)
968#endif
969 END SUBROUTINE gauxc_check_status
970
971!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
972! From here on, it's private helpers !
973!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
974
975#ifdef __GAUXC
976
977! **************************************************************************************************
978!> \brief ...
979!> \param spec ...
980!> \return ...
981! **************************************************************************************************
982 PURE FUNCTION read_execution_space(spec) RESULT(val)
983 CHARACTER(len=*), INTENT(IN) :: spec
984 INTEGER(c_int) :: val
985
986 CHARACTER(len=LEN(spec)) :: spec_upper
987
988 spec_upper = spec
989 CALL uppercase(spec_upper)
990
991 SELECT CASE (spec_upper)
992 CASE ("HOST")
993 val = gauxc_executionspace%host
994 CASE ("DEVICE")
995 val = gauxc_executionspace%device
996 CASE DEFAULT
997 val = gauxc_executionspace%host
998 END SELECT
999 END FUNCTION read_execution_space
1000
1001! **************************************************************************************************
1002!> \brief ...
1003!> \param spec ...
1004!> \return ...
1005! **************************************************************************************************
1006 PURE FUNCTION read_atomic_grid_size(spec) RESULT(val)
1007 CHARACTER(len=*), INTENT(IN) :: spec
1008 INTEGER(c_int) :: val
1009
1010 CHARACTER(len=LEN(spec)) :: spec_upper
1011
1012 spec_upper = spec
1013 CALL uppercase(spec_upper)
1014
1015 SELECT CASE (spec_upper)
1016 CASE ("FINE")
1017 val = gauxc_atomicgridsizedefault%finegrid
1018 CASE ("ULTRAFINE")
1019 val = gauxc_atomicgridsizedefault%ultrafinegrid
1020 CASE ("SUPERFINE")
1021 val = gauxc_atomicgridsizedefault%superfinegrid
1022 CASE ("GM3")
1023 val = gauxc_atomicgridsizedefault%gm3
1024 CASE ("GM5")
1025 val = gauxc_atomicgridsizedefault%gm5
1026 CASE DEFAULT
1027 val = gauxc_atomicgridsizedefault%finegrid
1028 END SELECT
1029 END FUNCTION read_atomic_grid_size
1030
1031! **************************************************************************************************
1032!> \brief ...
1033!> \param spec ...
1034!> \return ...
1035! **************************************************************************************************
1036 PURE FUNCTION read_radial_quad(spec) RESULT(val)
1037 CHARACTER(len=*), INTENT(IN) :: spec
1038 INTEGER(c_int) :: val
1039
1040 CHARACTER(len=LEN(spec)) :: spec_upper
1041
1042 spec_upper = spec
1043 CALL uppercase(spec_upper)
1044
1045 SELECT CASE (spec_upper)
1046 CASE ("BECKE")
1047 val = gauxc_radialquad%becke
1048 CASE ("MURAKNOWLES")
1049 val = gauxc_radialquad%mura_knowles
1050 CASE ("TREUTLERAHLRICHS")
1051 val = gauxc_radialquad%treutler_ahlrichs
1052 CASE ("MURRAYHANDYLAMING")
1053 val = gauxc_radialquad%murray_handy_laming
1054 CASE DEFAULT
1055 val = gauxc_radialquad%mura_knowles
1056 END SELECT
1057 END FUNCTION read_radial_quad
1058
1059! **************************************************************************************************
1060!> \brief ...
1061!> \param spec ...
1062!> \return ...
1063! **************************************************************************************************
1064 PURE FUNCTION read_pruning_scheme(spec) RESULT(val)
1065 CHARACTER(len=*), INTENT(IN) :: spec
1066 INTEGER(c_int) :: val
1067
1068 CHARACTER(len=LEN(spec)) :: spec_upper
1069
1070 spec_upper = spec
1071 CALL uppercase(spec_upper)
1072
1073 SELECT CASE (spec_upper)
1074 CASE ("UNPRUNED")
1075 val = gauxc_pruningscheme%unpruned
1076 CASE ("ROBUST")
1077 val = gauxc_pruningscheme%robust
1078 CASE ("TREUTLER")
1079 val = gauxc_pruningscheme%treutler
1080 CASE DEFAULT
1081 val = gauxc_pruningscheme%robust
1082 END SELECT
1083 END FUNCTION read_pruning_scheme
1084
1085#endif
1086
1087END MODULE xc_gauxc_interface
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::write_gto_basis_set
subroutine, public write_gto_basis_set(gto_basis_set, output_unit, header)
Write a Gaussian-type orbital (GTO) basis set data set to the output unit.
Definition basis_set_types.F:1867
cp_dbcsr_api
Definition cp_dbcsr_api.F:8
integrator
Provides integrator routines (velocity verlet) for all the ensemble types.
Definition integrator.F:26
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::get_ptable_info
subroutine, public get_ptable_info(symbol, number, amass, ielement, covalent_radius, metallic_radius, vdw_radius, found)
Pass information about the kind given the element symbol.
Definition periodic_table.F:85
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition physcon.F:147
qs_integral_utils
Some utility functions for the calculation of integrals.
Definition qs_integral_utils.F:14
qs_integral_utils::basis_set_list_setup
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Definition qs_integral_utils.F:149
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:468
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_dbcsr_api::dbcsr_p_type
Definition cp_dbcsr_api.F:172
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:177
print_gauxc_status_message
subroutine print_gauxc_status_message(status)
...
Definition xc_gauxc_interface.F:226