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 : ! **************************************************************************************************
9 : !> \brief Contains utility routines for the active space module
10 : !> \par History
11 : !> 04.2023 created [SB]
12 : !> \author SB
13 : ! **************************************************************************************************
14 : MODULE qs_active_space_utils
15 :
16 : USE cp_dbcsr_api, ONLY: dbcsr_csr_type
17 : USE cp_fm_types, ONLY: cp_fm_get_element,&
18 : cp_fm_get_info,&
19 : cp_fm_type
20 : USE kinds, ONLY: dp
21 : USE message_passing, ONLY: mp_comm_type
22 : USE qs_active_space_types, ONLY: csr_idx_from_combined,&
23 : csr_idx_to_combined,&
24 : eri_type
25 : #include "./base/base_uses.f90"
26 :
27 : IMPLICIT NONE
28 :
29 : PRIVATE
30 :
31 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_active_space_utils'
32 :
33 : PUBLIC :: subspace_matrix_to_array, eri_to_array, find_degenerate_splits
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief Copy a (square portion) of a `cp_fm_type` matrix to a standard 1D Fortran array
39 : !> \param source_matrix the matrix from where the data is taken
40 : !> \param target_array the array were the data is copied to
41 : !> \param row_index a list containing the row subspace indices
42 : !> \param col_index a list containing the column subspace indices
43 : ! **************************************************************************************************
44 54 : SUBROUTINE subspace_matrix_to_array(source_matrix, target_array, row_index, col_index)
45 : TYPE(cp_fm_type), INTENT(IN) :: source_matrix
46 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: target_array
47 : INTEGER, DIMENSION(:), INTENT(IN) :: row_index, col_index
48 :
49 : INTEGER :: i, i_sub, j, j_sub, max_col, max_row, &
50 : ncols, nrows
51 : REAL(KIND=dp) :: mval
52 :
53 54 : CALL cp_fm_get_info(source_matrix, nrow_global=max_row, ncol_global=max_col)
54 54 : nrows = SIZE(row_index)
55 54 : ncols = SIZE(col_index)
56 :
57 304 : CPASSERT(MAXVAL(row_index) <= max_row)
58 280 : CPASSERT(MAXVAL(col_index) <= max_col)
59 304 : CPASSERT(MINVAL(row_index) > 0)
60 280 : CPASSERT(MINVAL(col_index) > 0)
61 54 : CPASSERT(nrows <= max_row)
62 54 : CPASSERT(ncols <= max_col)
63 :
64 54 : CPASSERT(SIZE(target_array) == nrows*ncols)
65 :
66 280 : DO j = 1, ncols
67 226 : j_sub = col_index(j)
68 1446 : DO i = 1, nrows
69 1166 : i_sub = row_index(i)
70 1166 : CALL cp_fm_get_element(source_matrix, i_sub, j_sub, mval)
71 1392 : target_array(i + (j - 1)*nrows) = mval
72 : END DO
73 : END DO
74 54 : END SUBROUTINE subspace_matrix_to_array
75 :
76 : ! **************************************************************************************************
77 : !> \brief Copy the eri tensor for spins isp1 and isp2 to a standard 1D Fortran array
78 : !> \param eri_env the eri environment
79 : !> \param array the 1D Fortran array where the eri are copied to
80 : !> \param active_orbitals a list containing the active orbitals indices
81 : !> \param spin1 the spin of the bra
82 : !> \param spin2 the spin of the ket
83 : ! **************************************************************************************************
84 40 : SUBROUTINE eri_to_array(eri_env, array, active_orbitals, spin1, spin2)
85 : TYPE(eri_type), INTENT(IN) :: eri_env
86 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: array
87 : INTEGER, DIMENSION(:, :), INTENT(IN) :: active_orbitals
88 : INTEGER, INTENT(IN) :: spin1, spin2
89 :
90 : INTEGER :: i, i1, i12, i12l, i2, i3, i34, i34l, i4, &
91 : ijkl, ijlk, irptr, j, jikl, jilk, k, &
92 : klij, klji, l, lkij, lkji, nindex, &
93 : nmo_active, nmo_max
94 : REAL(KIND=dp) :: erival
95 : TYPE(dbcsr_csr_type), POINTER :: eri
96 : TYPE(mp_comm_type) :: mp_group
97 :
98 40 : nmo_active = SIZE(active_orbitals, 1)
99 40 : nmo_max = eri_env%norb
100 40 : nindex = (nmo_max*(nmo_max + 1))/2
101 34 : IF (spin1 == 1 .AND. spin2 == 1) THEN
102 28 : eri => eri_env%eri(1)%csr_mat
103 12 : ELSE IF ((spin1 == 1 .AND. spin2 == 2) .OR. (spin1 == 2 .AND. spin2 == 1)) THEN
104 6 : eri => eri_env%eri(2)%csr_mat
105 : ELSE
106 6 : eri => eri_env%eri(3)%csr_mat
107 : END IF
108 :
109 40 : CALL mp_group%set_handle(eri%mp_group%get_handle())
110 :
111 19850 : array = 0.0_dp
112 :
113 214 : DO i = 1, nmo_active
114 174 : i1 = active_orbitals(i, spin1)
115 702 : DO j = i, nmo_active
116 488 : i2 = active_orbitals(j, spin1)
117 488 : i12 = csr_idx_to_combined(i1, i2, nmo_max)
118 662 : IF (MOD(i12 - 1, eri_env%comm_exchange%num_pe) == eri_env%comm_exchange%mepos) THEN
119 488 : i12l = (i12 - 1)/eri_env%comm_exchange%num_pe + 1
120 488 : irptr = eri%rowptr_local(i12l) - 1
121 2131 : DO i34l = 1, eri%nzerow_local(i12l)
122 1643 : i34 = eri%colind_local(irptr + i34l)
123 1643 : CALL csr_idx_from_combined(i34, nmo_max, i3, i4)
124 4755 : k = FINDLOC(active_orbitals(:, spin2), i3, dim=1)
125 6568 : l = FINDLOC(active_orbitals(:, spin2), i4, dim=1)
126 1643 : erival = eri%nzval_local%r_dp(irptr + i34l)
127 :
128 : ! 8-fold permutational symmetry
129 1643 : ijkl = i + (j - 1)*nmo_active + (k - 1)*nmo_active**2 + (l - 1)*nmo_active**3
130 1643 : jikl = j + (i - 1)*nmo_active + (k - 1)*nmo_active**2 + (l - 1)*nmo_active**3
131 1643 : ijlk = i + (j - 1)*nmo_active + (l - 1)*nmo_active**2 + (k - 1)*nmo_active**3
132 1643 : jilk = j + (i - 1)*nmo_active + (l - 1)*nmo_active**2 + (k - 1)*nmo_active**3
133 1643 : array(ijkl) = erival
134 1643 : array(jikl) = erival
135 1643 : array(ijlk) = erival
136 1643 : array(jilk) = erival
137 3774 : IF (spin1 == spin2) THEN
138 1343 : klij = k + (l - 1)*nmo_active + (i - 1)*nmo_active**2 + (j - 1)*nmo_active**3
139 1343 : lkij = l + (k - 1)*nmo_active + (i - 1)*nmo_active**2 + (j - 1)*nmo_active**3
140 1343 : klji = k + (l - 1)*nmo_active + (j - 1)*nmo_active**2 + (i - 1)*nmo_active**3
141 1343 : lkji = l + (k - 1)*nmo_active + (j - 1)*nmo_active**2 + (i - 1)*nmo_active**3
142 1343 : array(klij) = erival
143 1343 : array(lkij) = erival
144 1343 : array(klji) = erival
145 1343 : array(lkji) = erival
146 : END IF
147 : END DO
148 : END IF
149 : END DO
150 : END DO
151 39660 : CALL mp_group%sum(array)
152 :
153 40 : END SUBROUTINE eri_to_array
154 :
155 : ! **************************************************************************************************
156 : !> \brief Finds pairs of degenerate orbitals of which one is active and the other is not, i.e.
157 : !> sets of degenerate orbitals that are split by the active space.
158 : !> \param eigenvalues orbital energies of the available orbitals
159 : !> \param active indices of the active orbitals; indices outside the available orbitals are skipped
160 : !> \param eps_degenerate energy difference below which two orbitals count as degenerate
161 : !> \param npairs number of split pairs
162 : !> \param first_active active orbital of the first split pair (0 if there is none)
163 : !> \param first_other non-active orbital of the first split pair (0 if there is none)
164 : ! **************************************************************************************************
165 124 : PURE SUBROUTINE find_degenerate_splits(eigenvalues, active, eps_degenerate, npairs, &
166 : first_active, first_other)
167 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
168 : INTEGER, DIMENSION(:), INTENT(IN) :: active
169 : REAL(KIND=dp), INTENT(IN) :: eps_degenerate
170 : INTEGER, INTENT(OUT) :: npairs, first_active, first_other
171 :
172 : INTEGER :: i, k, m
173 :
174 124 : npairs = 0
175 124 : first_active = 0
176 124 : first_other = 0
177 498 : DO i = 1, SIZE(active)
178 374 : m = active(i)
179 374 : IF (m < 1 .OR. m > SIZE(eigenvalues)) CYCLE
180 2946 : DO k = 1, SIZE(eigenvalues)
181 8332 : IF (ANY(active == k)) CYCLE
182 1592 : IF (ABS(eigenvalues(k) - eigenvalues(m)) < eps_degenerate) THEN
183 0 : npairs = npairs + 1
184 0 : IF (npairs == 1) THEN
185 0 : first_active = m
186 0 : first_other = k
187 : END IF
188 : END IF
189 : END DO
190 : END DO
191 :
192 124 : END SUBROUTINE find_degenerate_splits
193 :
194 : END MODULE qs_active_space_utils
|
