Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief The types needed for the calculation of active space Hamiltonians
10 : !> \par History
11 : !> 04.2016 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_active_space_types
15 :
16 : USE cp_dbcsr_api, ONLY: dbcsr_csr_destroy,&
17 : dbcsr_csr_p_type,&
18 : dbcsr_p_type
19 : USE cp_dbcsr_operations, ONLY: dbcsr_deallocate_matrix_set
20 : USE cp_fm_types, ONLY: cp_fm_release,&
21 : cp_fm_type
22 : USE input_constants, ONLY: eri_method_gpw_ht
23 : USE kinds, ONLY: default_path_length,&
24 : dp
25 : USE message_passing, ONLY: mp_comm_type
26 : USE qs_mo_types, ONLY: deallocate_mo_set,&
27 : mo_set_type
28 : #include "./base/base_uses.f90"
29 :
30 : IMPLICIT NONE
31 : PRIVATE
32 :
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_active_space_types'
34 :
35 : PUBLIC :: active_space_type, eri_type, eri_type_eri_element_func
36 : PUBLIC :: create_active_space_type, release_active_space_type
37 : PUBLIC :: csr_idx_to_combined, csr_idx_from_combined, get_irange_csr
38 :
39 : ! **************************************************************************************************
40 : !> \brief Quantities needed for AS determination
41 : !> \author JGH
42 : ! **************************************************************************************************
43 : TYPE eri_gpw_type
44 : LOGICAL :: redo_poisson = .FALSE.
45 : LOGICAL :: store_wfn = .FALSE.
46 : REAL(KIND=dp) :: cutoff = 0.0_dp
47 : REAL(KIND=dp) :: rel_cutoff = 0.0_dp
48 : REAL(KIND=dp) :: eps_grid = 0.0_dp
49 : REAL(KIND=dp) :: eps_filter = 0.0_dp
50 : INTEGER :: print_level = 0
51 : INTEGER :: group_size = 0
52 : END TYPE eri_gpw_type
53 :
54 : TYPE eri_type
55 : INTEGER :: method = 0
56 : INTEGER :: OPERATOR = 0
57 : REAL(KIND=dp) :: operator_parameter = 0.0_dp
58 : INTEGER, DIMENSION(3) :: periodicity = 0
59 : REAL(KIND=dp) :: cutoff_radius = 0.0_dp
60 : REAL(KIND=dp) :: eps_integral = 0.0_dp
61 : TYPE(eri_gpw_type) :: eri_gpw = eri_gpw_type()
62 : TYPE(dbcsr_csr_p_type), &
63 : DIMENSION(:), POINTER :: eri => NULL()
64 : INTEGER :: norb = 0
65 :
66 : CONTAINS
67 : PROCEDURE :: eri_foreach => eri_type_eri_foreach
68 : END TYPE eri_type
69 :
70 : ! **************************************************************************************************
71 : !> \brief Abstract function object for the `eri_type_eri_foreach` method
72 : ! **************************************************************************************************
73 : TYPE, ABSTRACT :: eri_type_eri_element_func
74 : CONTAINS
75 : PROCEDURE(eri_type_eri_element_func_interface), DEFERRED :: func
76 : END TYPE eri_type_eri_element_func
77 :
78 : TYPE active_space_type
79 : INTEGER :: nelec_active = 0
80 : INTEGER :: nelec_inactive = 0
81 : INTEGER :: nelec_total = 0
82 : INTEGER, POINTER, DIMENSION(:, :) :: active_orbitals => NULL()
83 : INTEGER, POINTER, DIMENSION(:, :) :: inactive_orbitals => NULL()
84 : INTEGER :: nmo_active = 0
85 : INTEGER :: nmo_inactive = 0
86 : INTEGER :: multiplicity = 0
87 : INTEGER :: nspins = 0
88 : LOGICAL :: molecule = .FALSE.
89 : INTEGER :: model = 0
90 : REAL(KIND=dp) :: energy_total = 0.0_dp
91 : REAL(KIND=dp) :: energy_ref = 0.0_dp
92 : REAL(KIND=dp) :: energy_inactive = 0.0_dp
93 : REAL(KIND=dp) :: energy_active = 0.0_dp
94 : LOGICAL :: do_scf_embedding = .FALSE.
95 : LOGICAL :: qcschema = .FALSE.
96 : LOGICAL :: fcidump = .FALSE.
97 : CHARACTER(LEN=default_path_length) :: qcschema_filename = ''
98 : TYPE(eri_type) :: eri = eri_type()
99 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_active => NULL()
100 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos_inactive => NULL()
101 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: p_active => NULL()
102 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: ks_sub => NULL()
103 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: vxc_sub => NULL()
104 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: h_sub => NULL()
105 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: fock_sub => NULL()
106 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: pmat_inactive => NULL()
107 : END TYPE active_space_type
108 :
109 : ABSTRACT INTERFACE
110 : ! **************************************************************************************************
111 : !> \brief The function signature to be implemented by a child of `eri_type_eri_element_func`
112 : !> \param this object reference
113 : !> \param i i-index
114 : !> \param j j-index
115 : !> \param k k-index
116 : !> \param l l-index
117 : !> \param val value of the integral at (i,j,k.l)
118 : !> \return True if the ERI foreach loop should continue, false, if not
119 : ! **************************************************************************************************
120 : LOGICAL FUNCTION eri_type_eri_element_func_interface(this, i, j, k, l, val)
121 : IMPORT :: eri_type_eri_element_func, dp
122 : CLASS(eri_type_eri_element_func), INTENT(inout) :: this
123 : INTEGER, INTENT(in) :: i, j, k, l
124 : REAL(KIND=dp), INTENT(in) :: val
125 : END FUNCTION eri_type_eri_element_func_interface
126 : END INTERFACE
127 :
128 : ! **************************************************************************************************
129 :
130 : CONTAINS
131 :
132 : ! **************************************************************************************************
133 : !> \brief Creates an active space environment type, nullifying all quantities.
134 : !> \param active_space_env the active space environment to be initialized
135 : ! **************************************************************************************************
136 62 : SUBROUTINE create_active_space_type(active_space_env)
137 : TYPE(active_space_type), POINTER :: active_space_env
138 :
139 62 : IF (ASSOCIATED(active_space_env)) THEN
140 0 : CALL release_active_space_type(active_space_env)
141 : END IF
142 :
143 248 : ALLOCATE (active_space_env)
144 : NULLIFY (active_space_env%active_orbitals, active_space_env%inactive_orbitals)
145 : NULLIFY (active_space_env%mos_active, active_space_env%mos_inactive)
146 : NULLIFY (active_space_env%ks_sub, active_space_env%p_active)
147 : NULLIFY (active_space_env%vxc_sub, active_space_env%h_sub)
148 : NULLIFY (active_space_env%fock_sub, active_space_env%pmat_inactive)
149 :
150 62 : END SUBROUTINE create_active_space_type
151 :
152 : ! **************************************************************************************************
153 : !> \brief Releases all quantities in the active space environment.
154 : !> \param active_space_env the active space environment to be released
155 : ! **************************************************************************************************
156 62 : SUBROUTINE release_active_space_type(active_space_env)
157 : TYPE(active_space_type), POINTER :: active_space_env
158 :
159 : INTEGER :: imo
160 :
161 62 : IF (ASSOCIATED(active_space_env)) THEN
162 :
163 62 : IF (ASSOCIATED(active_space_env%active_orbitals)) THEN
164 62 : DEALLOCATE (active_space_env%active_orbitals)
165 : END IF
166 :
167 62 : IF (ASSOCIATED(active_space_env%inactive_orbitals)) THEN
168 62 : DEALLOCATE (active_space_env%inactive_orbitals)
169 : END IF
170 :
171 62 : IF (ASSOCIATED(active_space_env%mos_active)) THEN
172 136 : DO imo = 1, SIZE(active_space_env%mos_active)
173 136 : CALL deallocate_mo_set(active_space_env%mos_active(imo))
174 : END DO
175 62 : DEALLOCATE (active_space_env%mos_active)
176 : END IF
177 :
178 62 : IF (ASSOCIATED(active_space_env%mos_inactive)) THEN
179 136 : DO imo = 1, SIZE(active_space_env%mos_inactive)
180 136 : CALL deallocate_mo_set(active_space_env%mos_inactive(imo))
181 : END DO
182 62 : DEALLOCATE (active_space_env%mos_inactive)
183 : END IF
184 :
185 62 : CALL release_eri_type(active_space_env%eri)
186 :
187 62 : CALL cp_fm_release(active_space_env%p_active)
188 62 : CALL cp_fm_release(active_space_env%ks_sub)
189 62 : CALL cp_fm_release(active_space_env%vxc_sub)
190 62 : CALL cp_fm_release(active_space_env%h_sub)
191 62 : CALL cp_fm_release(active_space_env%fock_sub)
192 :
193 62 : IF (ASSOCIATED(active_space_env%pmat_inactive)) &
194 62 : CALL dbcsr_deallocate_matrix_set(active_space_env%pmat_inactive)
195 :
196 62 : DEALLOCATE (active_space_env)
197 : END IF
198 :
199 62 : END SUBROUTINE release_active_space_type
200 :
201 : ! **************************************************************************************************
202 : !> \brief Releases the ERI environment type.
203 : !> \param eri_env the ERI environment to be released
204 : ! **************************************************************************************************
205 62 : SUBROUTINE release_eri_type(eri_env)
206 : TYPE(eri_type) :: eri_env
207 :
208 : INTEGER :: i
209 :
210 62 : IF (ASSOCIATED(eri_env%eri)) THEN
211 :
212 148 : DO i = 1, SIZE(eri_env%eri)
213 86 : CALL dbcsr_csr_destroy(eri_env%eri(i)%csr_mat)
214 148 : DEALLOCATE (eri_env%eri(i)%csr_mat)
215 : END DO
216 62 : DEALLOCATE (eri_env%eri)
217 :
218 : END IF
219 :
220 62 : END SUBROUTINE release_eri_type
221 :
222 : ! **************************************************************************************************
223 : !> \brief calculates combined index (ij)
224 : !> \param i Index j
225 : !> \param j Index i
226 : !> \param n Dimension in i or j direction
227 : !> \returns The combined index
228 : !> \par History
229 : !> 04.2016 created [JGH]
230 : ! **************************************************************************************************
231 3396 : INTEGER FUNCTION csr_idx_to_combined(i, j, n) RESULT(ij)
232 : INTEGER, INTENT(IN) :: i, j, n
233 :
234 3396 : CPASSERT(i <= j)
235 3396 : CPASSERT(i <= n)
236 3396 : CPASSERT(j <= n)
237 :
238 3396 : ij = (i - 1)*n - ((i - 1)*(i - 2))/2 + (j - i + 1)
239 :
240 3396 : CPASSERT(ij <= (n*(n + 1))/2 .AND. 0 <= ij)
241 :
242 3396 : END FUNCTION csr_idx_to_combined
243 :
244 : ! **************************************************************************************************
245 : !> \brief extracts indices i and j from combined index ij
246 : !> \param ij The combined index
247 : !> \param n Dimension in i or j direction
248 : !> \param i Resulting i index
249 : !> \param j Resulting j index
250 : !> \par History
251 : !> 04.2016 created [JGH]
252 : ! **************************************************************************************************
253 3669 : SUBROUTINE csr_idx_from_combined(ij, n, i, j)
254 : INTEGER, INTENT(IN) :: ij, n
255 : INTEGER, INTENT(OUT) :: i, j
256 :
257 : INTEGER :: m, m0
258 :
259 3669 : m = MAX(ij/n, 1)
260 9225 : DO i = m, n
261 9225 : m0 = (i - 1)*n - ((i - 1)*(i - 2))/2
262 9225 : j = ij - m0 + i - 1
263 9225 : IF (j <= n) EXIT
264 : END DO
265 :
266 3669 : CPASSERT(i > 0 .AND. i <= n)
267 3669 : CPASSERT(j > 0 .AND. j <= n)
268 3669 : CPASSERT(i <= j)
269 :
270 3669 : END SUBROUTINE csr_idx_from_combined
271 :
272 : ! **************************************************************************************************
273 : !> \brief calculates index range for processor in group mp_group
274 : !> \param nindex the number of indices
275 : !> \param mp_group message-passing group ID
276 : !> \return a range tuple
277 : !> \par History
278 : !> 04.2016 created [JGH]
279 : ! **************************************************************************************************
280 320 : FUNCTION get_irange_csr(nindex, mp_group) RESULT(irange)
281 : USE message_passing, ONLY: mp_comm_type
282 : INTEGER, INTENT(IN) :: nindex
283 :
284 : CLASS(mp_comm_type), INTENT(IN) :: mp_group
285 : INTEGER, DIMENSION(2) :: irange
286 :
287 : REAL(KIND=dp) :: rat
288 :
289 : ASSOCIATE (numtask => mp_group%num_pe, taskid => mp_group%mepos)
290 :
291 320 : IF (numtask == 1 .AND. taskid == 0) THEN
292 0 : irange(1) = 1
293 0 : irange(2) = nindex
294 320 : ELSEIF (numtask >= nindex) THEN
295 0 : IF (taskid >= nindex) THEN
296 0 : irange(1) = 1
297 0 : irange(2) = 0
298 : ELSE
299 0 : irange(1) = taskid + 1
300 0 : irange(2) = taskid + 1
301 : END IF
302 : ELSE
303 320 : rat = REAL(nindex, KIND=dp)/REAL(numtask, KIND=dp)
304 320 : irange(1) = NINT(rat*taskid) + 1
305 320 : irange(2) = NINT(rat*taskid + rat)
306 : END IF
307 : END ASSOCIATE
308 :
309 320 : END FUNCTION get_irange_csr
310 :
311 : ! **************************************************************************************************
312 : !> \brief Calls the provided function for each element in the ERI
313 : !> \param this object reference
314 : !> \param nspin The spin number
315 : !> \param active_orbitals the active orbital indices
316 : !> \param fobj The function object from which to call `func(i, j, k, l, val)`
317 : !> \param spin1 the first spin value
318 : !> \param spin2 the second spin value
319 : !> \par History
320 : !> 04.2016 created [JHU]
321 : !> 06.2016 factored out from qs_a_s_methods:fcidump [TMU]
322 : !> \note Calls MPI, must be executed on all ranks.
323 : ! **************************************************************************************************
324 172 : SUBROUTINE eri_type_eri_foreach(this, nspin, active_orbitals, fobj, spin1, spin2)
325 : CLASS(eri_type), INTENT(in) :: this
326 : CLASS(eri_type_eri_element_func) :: fobj
327 : INTEGER, DIMENSION(:, :), INTENT(IN) :: active_orbitals
328 : INTEGER, OPTIONAL :: spin1, spin2
329 : INTEGER :: i1, i12, i12l, i2, i3, i34, i34l, i4, m1, m2, m3, m4, &
330 : irange(2), irptr, nspin, nindex, nmo, proc, nonzero_elements_local
331 172 : INTEGER, ALLOCATABLE, DIMENSION(:) :: colind, offsets, nonzero_elements_global
332 172 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: erival
333 : REAL(KIND=dp) :: erint
334 : TYPE(mp_comm_type) :: mp_group
335 :
336 172 : IF (.NOT. PRESENT(spin1)) THEN
337 0 : spin1 = nspin
338 : END IF
339 172 : IF (.NOT. PRESENT(spin2)) THEN
340 0 : spin2 = nspin
341 : END IF
342 :
343 : ASSOCIATE (eri => this%eri(nspin)%csr_mat, norb => this%norb)
344 172 : nindex = (norb*(norb + 1))/2
345 172 : CALL mp_group%set_handle(eri%mp_group%get_handle())
346 172 : nmo = SIZE(active_orbitals, 1)
347 : ! Irrelevant in case of half-transformed integrals
348 172 : irange = get_irange_csr(nindex, mp_group)
349 860 : ALLOCATE (erival(nindex), colind(nindex))
350 :
351 172 : IF (this%method == eri_method_gpw_ht) THEN
352 : ALLOCATE (offsets(0:mp_group%num_pe - 1), &
353 96 : nonzero_elements_global(0:mp_group%num_pe - 1))
354 : END IF
355 :
356 780 : DO m1 = 1, nmo
357 436 : i1 = active_orbitals(m1, spin1)
358 1436 : DO m2 = m1, nmo
359 828 : i2 = active_orbitals(m2, spin1)
360 828 : i12 = csr_idx_to_combined(i1, i2, norb)
361 :
362 828 : IF (this%method == eri_method_gpw_ht) THEN
363 : ! In case of half-transformed integrals, every process might carry integrals of a row
364 : ! The number of integrals varies between processes and rows (related to the randomized
365 : ! distribution of matrix blocks)
366 :
367 : ! 1) Collect the amount of local data from each process
368 72 : nonzero_elements_local = eri%nzerow_local(i12)
369 72 : CALL mp_group%allgather(nonzero_elements_local, nonzero_elements_global)
370 :
371 : ! 2) Prepare arrays for communication (calculate the offsets and the total number of elements)
372 72 : offsets(0) = 0
373 144 : DO proc = 1, mp_group%num_pe - 1
374 144 : offsets(proc) = offsets(proc - 1) + nonzero_elements_global(proc - 1)
375 : END DO
376 72 : nindex = offsets(mp_group%num_pe - 1) + nonzero_elements_global(mp_group%num_pe - 1)
377 72 : irptr = eri%rowptr_local(i12)
378 :
379 : ! Exchange actual data
380 : CALL mp_group%allgatherv(eri%colind_local(irptr:irptr + nonzero_elements_local - 1), &
381 132 : colind(1:nindex), nonzero_elements_global, offsets)
382 : CALL mp_group%allgatherv(eri%nzval_local%r_dp(irptr:irptr + nonzero_elements_local - 1), &
383 132 : erival(1:nindex), nonzero_elements_global, offsets)
384 : ELSE
385 : ! Here, the rows are distributed among the processes such that each process
386 : ! carries all integral of a set of rows
387 756 : IF (i12 >= irange(1) .AND. i12 <= irange(2)) THEN
388 378 : i12l = i12 - irange(1) + 1
389 378 : irptr = eri%rowptr_local(i12l)
390 378 : nindex = eri%nzerow_local(i12l)
391 1738 : colind(1:nindex) = eri%colind_local(irptr:irptr + nindex - 1)
392 1738 : erival(1:nindex) = eri%nzval_local%r_dp(irptr:irptr + nindex - 1)
393 : ELSE
394 5250 : erival = 0.0_dp
395 5250 : colind = 0
396 378 : nindex = 0
397 : END IF
398 :
399 : ! Thus, a simple summation is sufficient
400 756 : CALL mp_group%sum(nindex)
401 756 : CALL mp_group%sum(colind(1:nindex))
402 756 : CALL mp_group%sum(erival(1:nindex))
403 : END IF
404 :
405 4104 : DO i34l = 1, nindex
406 2840 : i34 = colind(i34l)
407 2840 : erint = erival(i34l)
408 2840 : CALL csr_idx_from_combined(i34, norb, i3, i4)
409 :
410 6812 : DO m3 = 1, nmo
411 6812 : IF (active_orbitals(m3, spin2) == i3) THEN
412 : EXIT
413 : END IF
414 : END DO
415 :
416 9220 : DO m4 = 1, nmo
417 9220 : IF (active_orbitals(m4, spin2) == i4) THEN
418 : EXIT
419 : END IF
420 : END DO
421 :
422 : ! terminate the loop prematurely if the function returns false
423 6508 : IF (.NOT. fobj%func(m1, m2, m3, m4, erint)) RETURN
424 : END DO
425 :
426 : END DO
427 : END DO
428 : END ASSOCIATE
429 344 : END SUBROUTINE eri_type_eri_foreach
430 :
431 0 : END MODULE qs_active_space_types
|
