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