qs_active_space_types.F

Go to the documentation of this file.

!--------------------------------------------------------------------------------------------------!
!   CP2K: A general program to perform molecular dynamics simulations                              !
!   Copyright 2000-2026 CP2K developers group <https://cp2k.org>                                   !
!                                                                                                  !
!   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
!--------------------------------------------------------------------------------------------------!
 
! **************************************************************************************************
!> \brief The types needed for the calculation of active space Hamiltonians
!> \par History
!>      04.2016 created [JGH]
!> \author JGH
! **************************************************************************************************
MODULE qs_active_space_types
   USE cp_dbcsr_api,                    ONLY: dbcsr_csr_destroy,&
                                              dbcsr_csr_p_type,&
                                              dbcsr_p_type
   USE cp_dbcsr_operations,             ONLY: dbcsr_deallocate_matrix_set
   USE cp_fm_types,                     ONLY: cp_fm_release,&
                                              cp_fm_type
   USE input_section_types,             ONLY: section_vals_type
   USE kinds,                           ONLY: default_path_length,&
                                              dp
   USE message_passing,                 ONLY: mp_comm_null,&
                                              mp_comm_type,&
                                              mp_para_env_release,&
                                              mp_para_env_type
   USE qs_density_mixing_types,         ONLY: direct_mixing_nr,&
                                              mixing_storage_release,&
                                              mixing_storage_type
   USE qs_mo_types,                     ONLY: deallocate_mo_set,&
                                              mo_set_type
#include "./base/base_uses.f90"
 
   IMPLICIT NONE
   PRIVATE
 
   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_active_space_types'
 
   PUBLIC :: active_space_type, eri_type, eri_type_eri_element_func
   PUBLIC :: create_active_space_type, release_active_space_type
   PUBLIC :: csr_idx_to_combined, csr_idx_from_combined
 
! **************************************************************************************************
!> \brief Quantities needed for AS determination
!> \author JGH
! **************************************************************************************************
   TYPE eri_gpw_type
      LOGICAL                       :: redo_poisson = .false.
      LOGICAL                       :: store_wfn = .false.
      REAL(KIND=dp)                 :: cutoff = 0.0_dp
      REAL(KIND=dp)                 :: rel_cutoff = 0.0_dp
      REAL(KIND=dp)                 :: eps_grid = 0.0_dp
      REAL(KIND=dp)                 :: eps_filter = 0.0_dp
      INTEGER                       :: print_level = 0
      INTEGER                       :: group_size = 0
   END TYPE eri_gpw_type
 
   TYPE eri_type
      INTEGER                       :: method = 0
      INTEGER                       :: operator = 0
      LOGICAL                       :: enlarge_cell = .false.
      REAL(kind=dp)                 :: omega = 0.0_dp
      INTEGER, DIMENSION(3)         :: periodicity = 0
      REAL(kind=dp), DIMENSION(3)   :: eri_cell = 0
      REAL(kind=dp), DIMENSION(3)   :: eri_cell_angles = 0
      REAL(kind=dp)                 :: cutoff_radius = 0.0_dp
      REAL(kind=dp)                 :: eps_integral = 0.0_dp
      TYPE(eri_gpw_type)            :: eri_gpw = eri_gpw_type()
      TYPE(dbcsr_csr_p_type), &
         DIMENSION(:), POINTER      :: eri => null()
      INTEGER                       :: norb = 0
      TYPE(mp_para_env_type), POINTER :: para_env_sub => null()
      TYPE(mp_comm_type) :: comm_exchange = mp_comm_null
   CONTAINS
      PROCEDURE :: eri_foreach => eri_type_eri_foreach
   END TYPE eri_type
 
! **************************************************************************************************
!> \brief Abstract function object for the `eri_type_eri_foreach` method
! **************************************************************************************************
   TYPE, ABSTRACT :: eri_type_eri_element_func
   CONTAINS
      PROCEDURE(eri_type_eri_element_func_interface), DEFERRED :: func
   END TYPE eri_type_eri_element_func
 
   TYPE active_space_type
      INTEGER                                      :: nelec_active = 0
      INTEGER                                      :: nelec_inactive = 0
      INTEGER                                      :: nelec_total = 0
      INTEGER, POINTER, DIMENSION(:, :)            :: active_orbitals => null()
      INTEGER, POINTER, DIMENSION(:, :)            :: inactive_orbitals => null()
      INTEGER                                      :: nmo_active = 0
      INTEGER                                      :: nmo_inactive = 0
      INTEGER                                      :: multiplicity = 0
      INTEGER                                      :: nspins = 0
      LOGICAL                                      :: restricted_orbitals = .false.
      LOGICAL                                      :: molecule = .false.
      INTEGER                                      :: model = 0
      REAL(kind=dp)                                :: energy_total = 0.0_dp
      REAL(kind=dp)                                :: energy_ref = 0.0_dp
      REAL(kind=dp)                                :: energy_inactive = 0.0_dp
      REAL(kind=dp)                                :: energy_active = 0.0_dp
      REAL(kind=dp)                                :: eps_fci = 0.0_dp
      REAL(kind=dp)                                :: alpha = 0.0_dp
      INTEGER                                      :: as_mixing_method = direct_mixing_nr
      INTEGER                                      :: as_mixing_iter = 0
      INTEGER                                      :: as_mixing_dim = 0
      LOGICAL                                      :: do_scf_embedding = .false.
      LOGICAL                                      :: qcschema = .false.
      LOGICAL                                      :: fcidump = .false.
      CHARACTER(LEN=default_path_length)           :: qcschema_filename = ''
      TYPE(eri_type)                               :: eri = eri_type()
      REAL(kind=dp), DIMENSION(:), POINTER         :: as_mix_r_old => null()
      REAL(kind=dp), DIMENSION(:), POINTER         :: as_mix_weight => null()
      REAL(kind=dp), DIMENSION(:), POINTER         :: as_mix_x_old => null()
      REAL(kind=dp), DIMENSION(:, :), POINTER      :: as_mix_r_buffer => null()
      REAL(kind=dp), DIMENSION(:, :), POINTER      :: as_mix_x_buffer => null()
      TYPE(mixing_storage_type), POINTER           :: as_mixing_store => null()
      TYPE(mo_set_type), DIMENSION(:), POINTER     :: mos_active => null()
      TYPE(mo_set_type), DIMENSION(:), POINTER     :: mos_inactive => null()
      TYPE(cp_fm_type), DIMENSION(:), POINTER      :: p_active => null()
      TYPE(cp_fm_type), DIMENSION(:), POINTER      :: ks_sub => null()
      TYPE(cp_fm_type), DIMENSION(:), POINTER      :: vxc_sub => null()
      TYPE(cp_fm_type), DIMENSION(:), POINTER      :: h_sub => null()
      TYPE(cp_fm_type), DIMENSION(:), POINTER      :: fock_sub => null()
      TYPE(cp_fm_type), DIMENSION(:), POINTER      :: sab_sub => null()
      TYPE(dbcsr_p_type), DIMENSION(:), POINTER    :: pmat_inactive => null()
      TYPE(section_vals_type), POINTER             :: xc_section => null()
   END TYPE active_space_type
 
   abstract INTERFACE
! **************************************************************************************************
!> \brief The function signature to be implemented by a child of `eri_type_eri_element_func`
!> \param this object reference
!> \param i i-index
!> \param j j-index
!> \param k k-index
!> \param l l-index
!> \param val value of the integral at (i,j,k.l)
!> \return True if the ERI foreach loop should continue, false, if not
! **************************************************************************************************
      LOGICAL FUNCTION eri_type_eri_element_func_interface(this, i, j, k, l, val)
         IMPORT :: eri_type_eri_element_func, dp
         CLASS(eri_type_eri_element_func), INTENT(inout) :: this
         INTEGER, INTENT(in)                             :: i, j, k, l
         REAL(kind=dp), INTENT(in)                       :: val
      END FUNCTION eri_type_eri_element_func_interface
   END INTERFACE
 
! **************************************************************************************************
 
CONTAINS
 
! **************************************************************************************************
!> \brief Creates an active space environment type, nullifying all quantities.
!> \param active_space_env the active space environment to be initialized
! **************************************************************************************************
   SUBROUTINE create_active_space_type(active_space_env)
      TYPE(active_space_type), POINTER                   :: active_space_env
 
      IF (ASSOCIATED(active_space_env)) THEN
         CALL release_active_space_type(active_space_env)
      END IF
 
      ALLOCATE (active_space_env)
      NULLIFY (active_space_env%active_orbitals, active_space_env%inactive_orbitals)
      NULLIFY (active_space_env%mos_active, active_space_env%mos_inactive)
      NULLIFY (active_space_env%ks_sub, active_space_env%p_active)
      NULLIFY (active_space_env%vxc_sub, active_space_env%h_sub)
      NULLIFY (active_space_env%fock_sub, active_space_env%pmat_inactive)
      NULLIFY (active_space_env%as_mixing_store)
      NULLIFY (active_space_env%as_mix_r_old, active_space_env%as_mix_weight)
      NULLIFY (active_space_env%as_mix_x_old)
      NULLIFY (active_space_env%as_mix_r_buffer, active_space_env%as_mix_x_buffer)
 
   END SUBROUTINE create_active_space_type
 
! **************************************************************************************************
!> \brief Releases all quantities in the active space environment.
!> \param active_space_env the active space environment to be released
! **************************************************************************************************
   SUBROUTINE release_active_space_type(active_space_env)
      TYPE(active_space_type), POINTER                   :: active_space_env
 
      INTEGER                                            :: imo
 
      IF (ASSOCIATED(active_space_env)) THEN
 
         IF (ASSOCIATED(active_space_env%active_orbitals)) THEN
            DEALLOCATE (active_space_env%active_orbitals)
         END IF
 
         IF (ASSOCIATED(active_space_env%inactive_orbitals)) THEN
            DEALLOCATE (active_space_env%inactive_orbitals)
         END IF
 
         IF (ASSOCIATED(active_space_env%mos_active)) THEN
            DO imo = 1, SIZE(active_space_env%mos_active)
               CALL deallocate_mo_set(active_space_env%mos_active(imo))
            END DO
            DEALLOCATE (active_space_env%mos_active)
         END IF
 
         IF (ASSOCIATED(active_space_env%mos_inactive)) THEN
            DO imo = 1, SIZE(active_space_env%mos_inactive)
               CALL deallocate_mo_set(active_space_env%mos_inactive(imo))
            END DO
            DEALLOCATE (active_space_env%mos_inactive)
         END IF
 
         CALL release_eri_type(active_space_env%eri)
 
         CALL cp_fm_release(active_space_env%p_active)
         CALL cp_fm_release(active_space_env%ks_sub)
         CALL cp_fm_release(active_space_env%vxc_sub)
         CALL cp_fm_release(active_space_env%h_sub)
         CALL cp_fm_release(active_space_env%fock_sub)
         CALL cp_fm_release(active_space_env%sab_sub)
 
         IF (ASSOCIATED(active_space_env%as_mixing_store)) THEN
            CALL mixing_storage_release(active_space_env%as_mixing_store)
            DEALLOCATE (active_space_env%as_mixing_store)
         END IF
         IF (ASSOCIATED(active_space_env%as_mix_r_old)) THEN
            DEALLOCATE (active_space_env%as_mix_r_old)
         END IF
         IF (ASSOCIATED(active_space_env%as_mix_weight)) THEN
            DEALLOCATE (active_space_env%as_mix_weight)
         END IF
         IF (ASSOCIATED(active_space_env%as_mix_x_old)) THEN
            DEALLOCATE (active_space_env%as_mix_x_old)
         END IF
         IF (ASSOCIATED(active_space_env%as_mix_r_buffer)) THEN
            DEALLOCATE (active_space_env%as_mix_r_buffer)
         END IF
         IF (ASSOCIATED(active_space_env%as_mix_x_buffer)) THEN
            DEALLOCATE (active_space_env%as_mix_x_buffer)
         END IF
 
         IF (ASSOCIATED(active_space_env%pmat_inactive)) THEN
            CALL dbcsr_deallocate_matrix_set(active_space_env%pmat_inactive)
         END IF
 
         DEALLOCATE (active_space_env)
      END IF
 
   END SUBROUTINE release_active_space_type
 
! **************************************************************************************************
!> \brief Releases the ERI environment type.
!> \param eri_env the ERI environment to be released
! **************************************************************************************************
   SUBROUTINE release_eri_type(eri_env)
      TYPE(eri_type)                                     :: eri_env
 
      INTEGER                                            :: i
 
      IF (ASSOCIATED(eri_env%eri)) THEN
 
         DO i = 1, SIZE(eri_env%eri)
            CALL dbcsr_csr_destroy(eri_env%eri(i)%csr_mat)
            DEALLOCATE (eri_env%eri(i)%csr_mat)
         END DO
         CALL mp_para_env_release(eri_env%para_env_sub)
         CALL eri_env%comm_exchange%free()
         DEALLOCATE (eri_env%eri)
 
      END IF
 
   END SUBROUTINE release_eri_type
 
! **************************************************************************************************
!> \brief calculates combined index (ij)
!> \param i Index j
!> \param j Index i
!> \param n Dimension in i or j direction
!> \returns The combined index
!> \par History
!>      04.2016 created [JGH]
! **************************************************************************************************
   INTEGER FUNCTION csr_idx_to_combined(i, j, n) RESULT(ij)
      INTEGER, INTENT(IN)                                :: i, j, n
 
      cpassert(i <= j)
      cpassert(i <= n)
      cpassert(j <= n)
 
      ij = (i - 1)*n - ((i - 1)*(i - 2))/2 + (j - i + 1)
 
      cpassert(ij <= (n*(n + 1))/2 .AND. 0 <= ij)
 
   END FUNCTION csr_idx_to_combined
 
! **************************************************************************************************
!> \brief extracts indices i and j from combined index ij
!> \param ij The combined index
!> \param n Dimension in i or j direction
!> \param i Resulting i index
!> \param j Resulting j index
!> \par History
!>      04.2016 created [JGH]
! **************************************************************************************************
   SUBROUTINE csr_idx_from_combined(ij, n, i, j)
      INTEGER, INTENT(IN)                                :: ij, n
      INTEGER, INTENT(OUT)                               :: i, j
 
      INTEGER                                            :: m, m0
 
      m = max(ij/n, 1)
      DO i = m, n
         m0 = (i - 1)*n - ((i - 1)*(i - 2))/2
         j = ij - m0 + i - 1
         IF (j <= n) EXIT
      END DO
 
      cpassert(i > 0 .AND. i <= n)
      cpassert(j > 0 .AND. j <= n)
      cpassert(i <= j)
 
   END SUBROUTINE csr_idx_from_combined
 
! **************************************************************************************************
!> \brief Calls the provided function for each element in the ERI
!> \param this object reference
!> \param nspin The spin number
!> \param active_orbitals the active orbital indices
!> \param fobj The function object from which to call `func(i, j, k, l, val)`
!> \param spin1 the first spin value
!> \param spin2 the second spin value
!> \par History
!>      04.2016 created [JHU]
!>      06.2016 factored out from qs_a_s_methods:fcidump [TMU]
!> \note Calls MPI, must be executed on all ranks.
! **************************************************************************************************
   SUBROUTINE eri_type_eri_foreach(this, nspin, active_orbitals, fobj, spin1, spin2)
      CLASS(eri_type), INTENT(in)              :: this
      CLASS(eri_type_eri_element_func)         :: fobj
      INTEGER, DIMENSION(:, :), INTENT(IN)     :: active_orbitals
      INTEGER, OPTIONAL                        :: spin1, spin2
 
      CHARACTER(LEN=*), PARAMETER :: routinen = "eri_type_eri_foreach"
 
      INTEGER                                  :: i1, i12, i12l, i2, i3, i34, i34l, i4, m1, m2, m3, m4, &
                                                  irptr, nspin, nindex, nmo, proc, nonzero_elements_local, handle, dummy_int(1)
      INTEGER, ALLOCATABLE, DIMENSION(:)       :: colind, offsets, nonzero_elements_global
      REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: erival
      REAL(kind=dp)                            :: erint, dummy_real(1)
      TYPE(mp_comm_type) :: mp_group
 
      CALL timeset(routinen, handle)
 
      IF (.NOT. PRESENT(spin1)) THEN
         spin1 = nspin
      END IF
      IF (.NOT. PRESENT(spin2)) THEN
         spin2 = nspin
      END IF
 
      dummy_int = 0
      dummy_real = 0.0_dp
 
      associate(eri => this%eri(nspin)%csr_mat, norb => this%norb)
         nindex = (norb*(norb + 1))/2
         CALL mp_group%set_handle(eri%mp_group%get_handle())
         nmo = SIZE(active_orbitals, 1)
         ! Irrelevant in case of half-transformed integrals
         ALLOCATE (erival(nindex), colind(nindex))
         ALLOCATE (offsets(0:mp_group%num_pe - 1), &
                   nonzero_elements_global(0:mp_group%num_pe - 1))
 
         DO m1 = 1, nmo
            i1 = active_orbitals(m1, spin1)
            DO m2 = m1, nmo
               i2 = active_orbitals(m2, spin1)
               i12 = csr_idx_to_combined(i1, i2, norb)
               i12l = (i12 - 1)/this%comm_exchange%num_pe + 1
 
               ! In case of half-transformed integrals, every process might carry integrals of a row
               ! The number of integrals varies between processes and rows (related to the randomized
               ! distribution of matrix blocks)
 
               ! 1) Collect the amount of local data from each process
               nonzero_elements_local = 0
               IF (mod(i12 - 1, this%comm_exchange%num_pe) == this%comm_exchange%mepos) THEN
                  nonzero_elements_local = eri%nzerow_local(i12l)
               END IF
               CALL mp_group%allgather(nonzero_elements_local, nonzero_elements_global)
 
               ! 2) Prepare arrays for communication (calculate the offsets and the total number of elements)
               offsets(0) = 0
               DO proc = 1, mp_group%num_pe - 1
                  offsets(proc) = offsets(proc - 1) + nonzero_elements_global(proc - 1)
               END DO
               nindex = offsets(mp_group%num_pe - 1) + nonzero_elements_global(mp_group%num_pe - 1)
               irptr = 1
               IF (mod(i12 - 1, this%comm_exchange%num_pe) == this%comm_exchange%mepos) THEN
                  irptr = eri%rowptr_local(i12l)
 
                  ! Exchange actual data
                  CALL mp_group%allgatherv(eri%colind_local(irptr:irptr + nonzero_elements_local - 1), &
                                           colind(1:nindex), nonzero_elements_global, offsets)
                  CALL mp_group%allgatherv(eri%nzval_local%r_dp(irptr:irptr + nonzero_elements_local - 1), &
                                           erival(1:nindex), nonzero_elements_global, offsets)
               ELSE
                  CALL mp_group%allgatherv(dummy_int(1:0), colind(1:nindex), nonzero_elements_global, offsets)
                  CALL mp_group%allgatherv(dummy_real(1:0), erival(1:nindex), nonzero_elements_global, offsets)
               END IF
 
               DO i34l = 1, nindex
                  i34 = colind(i34l)
                  erint = erival(i34l)
                  CALL csr_idx_from_combined(i34, norb, i3, i4)
 
                  DO m3 = 1, nmo
                     IF (active_orbitals(m3, spin2) == i3) THEN
                        EXIT
                     END IF
                  END DO
 
                  DO m4 = 1, nmo
                     IF (active_orbitals(m4, spin2) == i4) THEN
                        EXIT
                     END IF
                  END DO
 
                  ! terminate the loop prematurely if the function returns false
                  IF (.NOT. fobj%func(m1, m2, m3, m4, erint)) RETURN
               END DO
 
            END DO
         END DO
      END associate
 
      CALL timestop(handle)
   END SUBROUTINE eri_type_eri_foreach
 
END MODULE qs_active_space_types