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