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