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 32 : 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 32 : CALL cp_fm_get_info(source_matrix, nrow_global=max_row, ncol_global=max_col)
54 32 : nrows = SIZE(row_index)
55 32 : ncols = SIZE(col_index)
56 :
57 150 : CPASSERT(MAXVAL(row_index) <= max_row)
58 126 : CPASSERT(MAXVAL(col_index) <= max_col)
59 150 : CPASSERT(MINVAL(row_index) > 0)
60 126 : CPASSERT(MINVAL(col_index) > 0)
61 32 : CPASSERT(nrows <= max_row)
62 32 : CPASSERT(ncols <= max_col)
63 :
64 32 : CPASSERT(SIZE(target_array) == nrows*ncols)
65 :
66 126 : DO j = 1, ncols
67 94 : j_sub = col_index(j)
68 612 : DO i = 1, nrows
69 486 : i_sub = row_index(i)
70 486 : CALL cp_fm_get_element(source_matrix, i_sub, j_sub, mval)
71 580 : target_array(i + (j - 1)*nrows) = mval
72 : END DO
73 : END DO
74 32 : 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 22 : 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 22 : nmo_active = SIZE(active_orbitals, 1)
99 22 : nmo_max = eri_env%norb
100 22 : nindex = (nmo_max*(nmo_max + 1))/2
101 18 : IF (spin1 == 1 .AND. spin2 == 1) THEN
102 14 : eri => eri_env%eri(1)%csr_mat
103 8 : ELSE IF ((spin1 == 1 .AND. spin2 == 2) .OR. (spin1 == 2 .AND. spin2 == 1)) THEN
104 4 : eri => eri_env%eri(2)%csr_mat
105 : ELSE
106 4 : eri => eri_env%eri(3)%csr_mat
107 : END IF
108 :
109 22 : CALL mp_group%set_handle(eri%mp_group%get_handle())
110 :
111 3064 : array = 0.0_dp
112 :
113 76 : DO i = 1, nmo_active
114 54 : i1 = active_orbitals(i, spin1)
115 184 : DO j = i, nmo_active
116 108 : i2 = active_orbitals(j, spin1)
117 108 : i12 = csr_idx_to_combined(i1, i2, nmo_max)
118 162 : IF (MOD(i12 - 1, eri_env%comm_exchange%num_pe) == eri_env%comm_exchange%mepos) THEN
119 108 : i12l = (i12 - 1)/eri_env%comm_exchange%num_pe + 1
120 108 : irptr = eri%rowptr_local(i12l) - 1
121 400 : DO i34l = 1, eri%nzerow_local(i12l)
122 292 : i34 = eri%colind_local(irptr + i34l)
123 292 : CALL csr_idx_from_combined(i34, nmo_max, i3, i4)
124 897 : k = FINDLOC(active_orbitals(:, spin2), i3, dim=1)
125 1248 : l = FINDLOC(active_orbitals(:, spin2), i4, dim=1)
126 292 : erival = eri%nzval_local%r_dp(irptr + i34l)
127 :
128 : ! 8-fold permutational symmetry
129 292 : ijkl = i + (j - 1)*nmo_active + (k - 1)*nmo_active**2 + (l - 1)*nmo_active**3
130 292 : jikl = j + (i - 1)*nmo_active + (k - 1)*nmo_active**2 + (l - 1)*nmo_active**3
131 292 : ijlk = i + (j - 1)*nmo_active + (l - 1)*nmo_active**2 + (k - 1)*nmo_active**3
132 292 : jilk = j + (i - 1)*nmo_active + (l - 1)*nmo_active**2 + (k - 1)*nmo_active**3
133 292 : array(ijkl) = erival
134 292 : array(jikl) = erival
135 292 : array(ijlk) = erival
136 292 : array(jilk) = erival
137 692 : IF (spin1 == spin2) THEN
138 282 : klij = k + (l - 1)*nmo_active + (i - 1)*nmo_active**2 + (j - 1)*nmo_active**3
139 282 : lkij = l + (k - 1)*nmo_active + (i - 1)*nmo_active**2 + (j - 1)*nmo_active**3
140 282 : klji = k + (l - 1)*nmo_active + (j - 1)*nmo_active**2 + (i - 1)*nmo_active**3
141 282 : lkji = l + (k - 1)*nmo_active + (j - 1)*nmo_active**2 + (i - 1)*nmo_active**3
142 282 : array(klij) = erival
143 282 : array(lkij) = erival
144 282 : array(klji) = erival
145 282 : array(lkji) = erival
146 : END IF
147 : END DO
148 : END IF
149 : END DO
150 : END DO
151 6106 : CALL mp_group%sum(array)
152 :
153 22 : 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 136 : 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 136 : npairs = 0
175 136 : first_active = 0
176 136 : first_other = 0
177 518 : DO i = 1, SIZE(active)
178 382 : m = active(i)
179 382 : IF (m < 1 .OR. m > SIZE(eigenvalues)) CYCLE
180 3010 : DO k = 1, SIZE(eigenvalues)
181 8420 : IF (ANY(active == k)) CYCLE
182 1704 : IF (ABS(eigenvalues(k) - eigenvalues(m)) < eps_degenerate) THEN
183 18 : npairs = npairs + 1
184 18 : IF (npairs == 1) THEN
185 6 : first_active = m
186 6 : first_other = k
187 : END IF
188 : END IF
189 : END DO
190 : END DO
191 :
192 136 : END SUBROUTINE find_degenerate_splits
193 :
194 : END MODULE qs_active_space_utils
|
