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 Interface between ALMO SCF and QS
10 : !> \par History
11 : !> 2011.05 created [Rustam Z Khaliullin]
12 : !> \author Rustam Z Khaliullin
13 : ! **************************************************************************************************
14 : MODULE almo_scf_qs
15 : USE almo_scf_types, ONLY: almo_mat_dim_aobasis,&
16 : almo_mat_dim_occ,&
17 : almo_mat_dim_virt,&
18 : almo_mat_dim_virt_disc,&
19 : almo_mat_dim_virt_full,&
20 : almo_scf_env_type
21 : USE atomic_kind_types, ONLY: get_atomic_kind
22 : USE cell_types, ONLY: cell_type,&
23 : pbc
24 : USE cp_control_types, ONLY: dft_control_type
25 : USE cp_dbcsr_api, ONLY: &
26 : dbcsr_complete_redistribute, dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, &
27 : dbcsr_distribution_get, dbcsr_distribution_new, dbcsr_distribution_release, &
28 : dbcsr_distribution_type, dbcsr_filter, dbcsr_finalize, dbcsr_get_block_p, dbcsr_get_info, &
29 : dbcsr_get_num_blocks, dbcsr_get_stored_coordinates, dbcsr_multiply, dbcsr_p_type, &
30 : dbcsr_put_block, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, &
31 : dbcsr_work_create
32 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
33 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
34 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
35 : cp_fm_struct_release,&
36 : cp_fm_struct_type
37 : USE cp_fm_types, ONLY: cp_fm_create,&
38 : cp_fm_release,&
39 : cp_fm_type
40 : USE cp_log_handling, ONLY: cp_get_default_logger,&
41 : cp_logger_get_default_unit_nr,&
42 : cp_logger_type
43 : USE cp_units, ONLY: cp_unit_to_cp2k
44 : USE input_constants, ONLY: almo_constraint_ao_overlap,&
45 : almo_constraint_block_diagonal,&
46 : almo_constraint_distance,&
47 : almo_domain_layout_molecular,&
48 : almo_mat_distr_atomic,&
49 : almo_mat_distr_molecular,&
50 : do_bondparm_covalent,&
51 : do_bondparm_vdw
52 : USE kinds, ONLY: dp
53 : USE message_passing, ONLY: mp_comm_type
54 : USE molecule_types, ONLY: get_molecule_set_info,&
55 : molecule_type
56 : USE particle_types, ONLY: particle_type
57 : USE qs_energy_types, ONLY: qs_energy_type
58 : USE qs_environment_types, ONLY: get_qs_env,&
59 : qs_environment_type,&
60 : set_qs_env
61 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
62 : USE qs_ks_types, ONLY: qs_ks_did_change,&
63 : qs_ks_env_type,&
64 : set_ks_env
65 : USE qs_mo_types, ONLY: allocate_mo_set,&
66 : deallocate_mo_set,&
67 : init_mo_set,&
68 : mo_set_type
69 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
70 : neighbor_list_iterate,&
71 : neighbor_list_iterator_create,&
72 : neighbor_list_iterator_p_type,&
73 : neighbor_list_iterator_release,&
74 : neighbor_list_set_p_type
75 : USE qs_rho_methods, ONLY: qs_rho_update_rho
76 : USE qs_rho_types, ONLY: qs_rho_get,&
77 : qs_rho_type
78 : USE qs_scf_types, ONLY: qs_scf_env_type,&
79 : scf_env_create
80 : #include "./base/base_uses.f90"
81 :
82 : IMPLICIT NONE
83 :
84 : PRIVATE
85 :
86 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'almo_scf_qs'
87 :
88 : PUBLIC :: matrix_almo_create, &
89 : almo_scf_construct_quencher, &
90 : calculate_w_matrix_almo, &
91 : init_almo_ks_matrix_via_qs, &
92 : almo_scf_update_ks_energy, &
93 : construct_qs_mos, &
94 : matrix_qs_to_almo, &
95 : almo_dm_to_almo_ks, &
96 : almo_dm_to_qs_env
97 :
98 : CONTAINS
99 :
100 : ! **************************************************************************************************
101 : !> \brief create the ALMO matrix templates
102 : !> \param matrix_new ...
103 : !> \param matrix_qs ...
104 : !> \param almo_scf_env ...
105 : !> \param name_new ...
106 : !> \param size_keys ...
107 : !> \param symmetry_new ...
108 : !> \param spin_key ...
109 : !> \param init_domains ...
110 : !> \par History
111 : !> 2011.05 created [Rustam Z Khaliullin]
112 : !> \author Rustam Z Khaliullin
113 : ! **************************************************************************************************
114 3596 : SUBROUTINE matrix_almo_create(matrix_new, matrix_qs, almo_scf_env, &
115 : name_new, size_keys, symmetry_new, &
116 : spin_key, init_domains)
117 :
118 : TYPE(dbcsr_type) :: matrix_new, matrix_qs
119 : TYPE(almo_scf_env_type), INTENT(IN) :: almo_scf_env
120 : CHARACTER(len=*), INTENT(IN) :: name_new
121 : INTEGER, DIMENSION(2), INTENT(IN) :: size_keys
122 : CHARACTER, INTENT(IN) :: symmetry_new
123 : INTEGER, INTENT(IN) :: spin_key
124 : LOGICAL, INTENT(IN) :: init_domains
125 :
126 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_almo_create'
127 :
128 : INTEGER :: dimen, handle, hold, iatom, iblock_col, &
129 : iblock_row, imol, mynode, natoms, &
130 : nblkrows_tot, nlength, nmols, row
131 3596 : INTEGER, DIMENSION(:), POINTER :: blk_distr, blk_sizes, block_sizes_new, col_blk_size, &
132 3596 : col_distr_new, col_sizes_new, distr_new_array, row_blk_size, row_distr_new, row_sizes_new
133 : LOGICAL :: active, one_dim_is_mo, tr
134 3596 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: new_block
135 : TYPE(dbcsr_distribution_type) :: dist_new, dist_qs
136 :
137 : ! dimension size: AO, MO, etc
138 : ! almo_mat_dim_aobasis - no. of AOs,
139 : ! almo_mat_dim_occ - no. of occupied MOs
140 : ! almo_mat_dim_domains - no. of domains
141 : ! symmetry type: dbcsr_type_no_symmetry, dbcsr_type_symmetric,
142 : ! dbcsr_type_antisymmetric, dbcsr_type_hermitian, dbcsr_type_antihermitian
143 : ! (see dbcsr_lib/dbcsr_types.F for other values)
144 : ! spin_key: either 1 or 2 (0 is allowed for matrics in the AO basis)
145 : ! TYPE(dbcsr_iterator_type) :: iter
146 : ! REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: allones
147 : !-----------------------------------------------------------------------
148 :
149 3596 : CALL timeset(routineN, handle)
150 :
151 : ! RZK-warning The structure of the matrices can be optimized:
152 : ! 1. Diagonal matrices must be distributed evenly over the processes.
153 : ! This can be achieved by distributing cpus: 012012-rows and 001122-cols
154 : ! block_diagonal_flag is introduced but not used
155 : ! 2. Multiplication of diagonally dominant matrices will be faster
156 : ! if the diagonal blocks are local to the same processes.
157 : ! 3. Systems of molecules of drastically different sizes might need
158 : ! better distribution.
159 :
160 : ! obtain distribution from the qs matrix - it might be useful
161 : ! to get the structure of the AO dimensions
162 3596 : CALL dbcsr_get_info(matrix_qs, distribution=dist_qs)
163 :
164 3596 : natoms = almo_scf_env%natoms
165 3596 : nmols = almo_scf_env%nmolecules
166 :
167 10788 : DO dimen = 1, 2 ! 1 - row, 2 - column dimension
168 :
169 : ! distribution pattern is the same for all matrix types (ao, occ, virt)
170 7192 : IF (dimen == 1) THEN !rows
171 3596 : CALL dbcsr_distribution_get(dist_qs, row_dist=blk_distr)
172 : ELSE !columns
173 3596 : CALL dbcsr_distribution_get(dist_qs, col_dist=blk_distr)
174 : END IF
175 :
176 7192 : IF (size_keys(dimen) == almo_mat_dim_aobasis) THEN ! this dimension is AO
177 :
178 : ! structure of an AO dimension can be copied from matrix_qs
179 2712 : CALL dbcsr_get_info(matrix_qs, row_blk_size=blk_sizes)
180 :
181 : ! atomic clustering of AOs
182 2712 : IF (almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN
183 0 : ALLOCATE (block_sizes_new(natoms), distr_new_array(natoms))
184 0 : block_sizes_new(:) = blk_sizes(:)
185 0 : distr_new_array(:) = blk_distr(:)
186 : ! molecular clustering of AOs
187 2712 : ELSE IF (almo_scf_env%mat_distr_aos == almo_mat_distr_molecular) THEN
188 10848 : ALLOCATE (block_sizes_new(nmols), distr_new_array(nmols))
189 20742 : block_sizes_new(:) = 0
190 47092 : DO iatom = 1, natoms
191 : block_sizes_new(almo_scf_env%domain_index_of_atom(iatom)) = &
192 : block_sizes_new(almo_scf_env%domain_index_of_atom(iatom)) + &
193 47092 : blk_sizes(iatom)
194 : END DO
195 20742 : DO imol = 1, nmols
196 : distr_new_array(imol) = &
197 20742 : blk_distr(almo_scf_env%first_atom_of_domain(imol))
198 : END DO
199 : ELSE
200 0 : CPABORT("Illegal distribution")
201 : END IF
202 :
203 : ELSE ! this dimension is not AO
204 :
205 : IF (size_keys(dimen) == almo_mat_dim_occ .OR. &
206 : size_keys(dimen) == almo_mat_dim_virt .OR. &
207 4480 : size_keys(dimen) == almo_mat_dim_virt_disc .OR. &
208 : size_keys(dimen) == almo_mat_dim_virt_full) THEN ! this dim is MO
209 :
210 : ! atomic clustering of MOs
211 4480 : IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN
212 0 : nlength = natoms
213 0 : ALLOCATE (block_sizes_new(nlength))
214 0 : block_sizes_new(:) = 0
215 : IF (size_keys(dimen) == almo_mat_dim_occ) THEN
216 : ! currently distributing atomic distr of mos is not allowed
217 : ! RZK-warning define nocc_of_atom and nvirt_atom to implement it
218 : !block_sizes_new(:)=almo_scf_env%nocc_of_atom(:,spin_key)
219 : ELSE IF (size_keys(dimen) == almo_mat_dim_virt) THEN
220 : !block_sizes_new(:)=almo_scf_env%nvirt_of_atom(:,spin_key)
221 : END IF
222 : ! molecular clustering of MOs
223 4480 : ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN
224 4480 : nlength = nmols
225 13440 : ALLOCATE (block_sizes_new(nlength))
226 4480 : IF (size_keys(dimen) == almo_mat_dim_occ) THEN
227 19240 : block_sizes_new(:) = almo_scf_env%nocc_of_domain(:, spin_key)
228 : ! Handle zero-electron fragments by adding one-orbital that
229 : ! must remain zero at all times
230 19240 : WHERE (block_sizes_new == 0) block_sizes_new = 1
231 1920 : ELSE IF (size_keys(dimen) == almo_mat_dim_virt_disc) THEN
232 0 : block_sizes_new(:) = almo_scf_env%nvirt_disc_of_domain(:, spin_key)
233 1920 : ELSE IF (size_keys(dimen) == almo_mat_dim_virt_full) THEN
234 5772 : block_sizes_new(:) = almo_scf_env%nvirt_full_of_domain(:, spin_key)
235 1152 : ELSE IF (size_keys(dimen) == almo_mat_dim_virt) THEN
236 8658 : block_sizes_new(:) = almo_scf_env%nvirt_of_domain(:, spin_key)
237 : END IF
238 : ELSE
239 0 : CPABORT("Illegal distribution")
240 : END IF
241 :
242 : ELSE
243 :
244 0 : CPABORT("Illegal dimension")
245 :
246 : END IF ! end choosing dim size (occ, virt)
247 :
248 : ! distribution for MOs is copied from AOs
249 13440 : ALLOCATE (distr_new_array(nlength))
250 : ! atomic clustering
251 4480 : IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN
252 0 : distr_new_array(:) = blk_distr(:)
253 : ! molecular clustering
254 4480 : ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN
255 33670 : DO imol = 1, nmols
256 : distr_new_array(imol) = &
257 33670 : blk_distr(almo_scf_env%first_atom_of_domain(imol))
258 : END DO
259 : END IF
260 : END IF ! end choosing dimension size (AOs vs .NOT.AOs)
261 :
262 : ! create final arrays
263 10788 : IF (dimen == 1) THEN !rows
264 3596 : row_sizes_new => block_sizes_new
265 3596 : row_distr_new => distr_new_array
266 : ELSE !columns
267 3596 : col_sizes_new => block_sizes_new
268 3596 : col_distr_new => distr_new_array
269 : END IF
270 : END DO ! both rows and columns are done
271 :
272 : ! Create the distribution
273 : CALL dbcsr_distribution_new(dist_new, template=dist_qs, &
274 : row_dist=row_distr_new, col_dist=col_distr_new, &
275 3596 : reuse_arrays=.TRUE.)
276 :
277 : ! Create the matrix
278 : CALL dbcsr_create(matrix_new, name_new, &
279 : dist_new, symmetry_new, &
280 3596 : row_sizes_new, col_sizes_new, reuse_arrays=.TRUE.)
281 3596 : CALL dbcsr_distribution_release(dist_new)
282 :
283 : ! fill out reqired blocks with 1.0_dp to tell the dbcsr library
284 : ! which blocks to keep
285 3596 : IF (init_domains) THEN
286 :
287 1426 : CALL dbcsr_distribution_get(dist_new, mynode=mynode)
288 1426 : CALL dbcsr_work_create(matrix_new, work_mutable=.TRUE.)
289 : CALL dbcsr_get_info(matrix_new, nblkrows_total=nblkrows_tot, &
290 1426 : row_blk_size=row_blk_size, col_blk_size=col_blk_size)
291 : ! start linear-scaling replacement:
292 : ! works only for molecular blocks AND molecular distributions
293 10870 : DO row = 1, nblkrows_tot
294 9444 : tr = .FALSE.
295 9444 : iblock_row = row
296 9444 : iblock_col = row
297 9444 : CALL dbcsr_get_stored_coordinates(matrix_new, iblock_row, iblock_col, hold)
298 :
299 10870 : IF (hold == mynode) THEN
300 :
301 14166 : active = .TRUE.
302 :
303 : one_dim_is_mo = .FALSE.
304 14166 : DO dimen = 1, 2 ! 1 - row, 2 - column dimension
305 14166 : IF (size_keys(dimen) == almo_mat_dim_occ) one_dim_is_mo = .TRUE.
306 : END DO
307 4722 : IF (one_dim_is_mo) THEN
308 1668 : IF (almo_scf_env%nocc_of_domain(row, spin_key) == 0) active = .FALSE.
309 : END IF
310 :
311 4722 : one_dim_is_mo = .FALSE.
312 14166 : DO dimen = 1, 2
313 14166 : IF (size_keys(dimen) == almo_mat_dim_virt) one_dim_is_mo = .TRUE.
314 : END DO
315 4722 : IF (one_dim_is_mo) THEN
316 417 : IF (almo_scf_env%nvirt_of_domain(row, spin_key) == 0) active = .FALSE.
317 : END IF
318 :
319 4722 : one_dim_is_mo = .FALSE.
320 14166 : DO dimen = 1, 2
321 14166 : IF (size_keys(dimen) == almo_mat_dim_virt_disc) one_dim_is_mo = .TRUE.
322 : END DO
323 4722 : IF (one_dim_is_mo) THEN
324 0 : IF (almo_scf_env%nvirt_disc_of_domain(row, spin_key) == 0) active = .FALSE.
325 : END IF
326 :
327 4722 : one_dim_is_mo = .FALSE.
328 14166 : DO dimen = 1, 2
329 14166 : IF (size_keys(dimen) == almo_mat_dim_virt_full) one_dim_is_mo = .TRUE.
330 : END DO
331 4722 : IF (one_dim_is_mo) THEN
332 417 : IF (almo_scf_env%nvirt_full_of_domain(row, spin_key) == 0) active = .FALSE.
333 : END IF
334 :
335 4688 : IF (active) THEN
336 17592 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
337 840230 : new_block(:, :) = 1.0_dp
338 4398 : CALL dbcsr_put_block(matrix_new, iblock_row, iblock_col, new_block)
339 4398 : DEALLOCATE (new_block)
340 : END IF
341 :
342 : END IF ! mynode
343 : END DO
344 : ! end lnear-scaling replacement
345 :
346 : END IF ! init_domains
347 :
348 3596 : CALL dbcsr_finalize(matrix_new)
349 :
350 3596 : CALL timestop(handle)
351 :
352 7192 : END SUBROUTINE matrix_almo_create
353 :
354 : ! **************************************************************************************************
355 : !> \brief convert between two types of matrices: QS style to ALMO style
356 : !> \param matrix_qs ...
357 : !> \param matrix_almo ...
358 : !> \param mat_distr_aos ...
359 : !> \par History
360 : !> 2011.06 created [Rustam Z Khaliullin]
361 : !> \author Rustam Z Khaliullin
362 : ! **************************************************************************************************
363 2080 : SUBROUTINE matrix_qs_to_almo(matrix_qs, matrix_almo, mat_distr_aos)
364 :
365 : TYPE(dbcsr_type) :: matrix_qs, matrix_almo
366 : INTEGER :: mat_distr_aos
367 :
368 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_qs_to_almo'
369 :
370 : INTEGER :: handle
371 : TYPE(dbcsr_type) :: matrix_qs_nosym
372 :
373 2080 : CALL timeset(routineN, handle)
374 : !RZK-warning if it's not a N(AO)xN(AO) matrix then stop
375 :
376 2080 : SELECT CASE (mat_distr_aos)
377 : CASE (almo_mat_distr_atomic)
378 : ! automatic data_type conversion
379 0 : CALL dbcsr_copy(matrix_almo, matrix_qs)
380 : CASE (almo_mat_distr_molecular)
381 : ! desymmetrize the qs matrix
382 2080 : CALL dbcsr_create(matrix_qs_nosym, template=matrix_qs, matrix_type=dbcsr_type_no_symmetry)
383 2080 : CALL dbcsr_desymmetrize(matrix_qs, matrix_qs_nosym)
384 :
385 : ! perform the magic complete_redistribute
386 : ! before calling complete_redistribute set all blocks to zero
387 : ! otherwise the non-zero elements of the redistributed matrix,
388 : ! which are in zero-blocks of the original matrix, will remain
389 : ! in the final redistributed matrix. this is a bug in
390 : ! complete_redistribute. RZK-warning it should be later corrected by calling
391 : ! dbcsr_set to 0.0 from within complete_redistribute
392 2080 : CALL dbcsr_set(matrix_almo, 0.0_dp)
393 2080 : CALL dbcsr_complete_redistribute(matrix_qs_nosym, matrix_almo)
394 2080 : CALL dbcsr_release(matrix_qs_nosym)
395 :
396 : CASE DEFAULT
397 2080 : CPABORT("Unknown mat_distr_aos for matrix_qs_to_almo")
398 : END SELECT
399 :
400 2080 : CALL timestop(handle)
401 :
402 2080 : END SUBROUTINE matrix_qs_to_almo
403 :
404 : ! **************************************************************************************************
405 : !> \brief convert between two types of matrices: ALMO style to QS style
406 : !> \param matrix_almo ...
407 : !> \param matrix_qs ...
408 : !> \param mat_distr_aos ...
409 : !> \par History
410 : !> 2011.06 created [Rustam Z Khaliullin]
411 : !> \author Rustam Z Khaliullin
412 : ! **************************************************************************************************
413 1862 : SUBROUTINE matrix_almo_to_qs(matrix_almo, matrix_qs, mat_distr_aos)
414 : TYPE(dbcsr_type) :: matrix_almo, matrix_qs
415 : INTEGER, INTENT(IN) :: mat_distr_aos
416 :
417 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_almo_to_qs'
418 :
419 : INTEGER :: handle
420 : TYPE(dbcsr_type) :: matrix_almo_redist
421 :
422 1862 : CALL timeset(routineN, handle)
423 : ! RZK-warning if it's not a N(AO)xN(AO) matrix then stop
424 :
425 1862 : SELECT CASE (mat_distr_aos)
426 : CASE (almo_mat_distr_atomic)
427 0 : CALL dbcsr_copy(matrix_qs, matrix_almo, keep_sparsity=.TRUE.)
428 : CASE (almo_mat_distr_molecular)
429 1862 : CALL dbcsr_create(matrix_almo_redist, template=matrix_qs)
430 1862 : CALL dbcsr_complete_redistribute(matrix_almo, matrix_almo_redist)
431 1862 : CALL dbcsr_set(matrix_qs, 0.0_dp)
432 1862 : CALL dbcsr_copy(matrix_qs, matrix_almo_redist, keep_sparsity=.TRUE.)
433 1862 : CALL dbcsr_release(matrix_almo_redist)
434 : CASE DEFAULT
435 1862 : CPABORT("Unknown mat_distr_aos for matrix_almo_to_qs")
436 : END SELECT
437 :
438 1862 : CALL timestop(handle)
439 :
440 1862 : END SUBROUTINE matrix_almo_to_qs
441 :
442 : ! **************************************************************************************************
443 : !> \brief Initialization of the QS and ALMO KS matrix
444 : !> \param qs_env ...
445 : !> \param matrix_ks ...
446 : !> \param mat_distr_aos ...
447 : !> \param eps_filter ...
448 : !> \par History
449 : !> 2011.05 created [Rustam Z Khaliullin]
450 : !> \author Rustam Z Khaliullin
451 : ! **************************************************************************************************
452 122 : SUBROUTINE init_almo_ks_matrix_via_qs(qs_env, matrix_ks, mat_distr_aos, eps_filter)
453 :
454 : TYPE(qs_environment_type), POINTER :: qs_env
455 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_ks
456 : INTEGER :: mat_distr_aos
457 : REAL(KIND=dp) :: eps_filter
458 :
459 : CHARACTER(len=*), PARAMETER :: routineN = 'init_almo_ks_matrix_via_qs'
460 :
461 : INTEGER :: handle, ispin, nspin
462 122 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_qs_ks, matrix_qs_s
463 : TYPE(dft_control_type), POINTER :: dft_control
464 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
465 122 : POINTER :: sab_orb
466 : TYPE(qs_ks_env_type), POINTER :: ks_env
467 :
468 122 : CALL timeset(routineN, handle)
469 :
470 122 : NULLIFY (sab_orb)
471 :
472 : ! get basic quantities from the qs_env
473 : CALL get_qs_env(qs_env, &
474 : dft_control=dft_control, &
475 : matrix_s=matrix_qs_s, &
476 : matrix_ks=matrix_qs_ks, &
477 : ks_env=ks_env, &
478 122 : sab_orb=sab_orb)
479 :
480 122 : nspin = dft_control%nspins
481 :
482 : ! create matrix_ks in the QS env if necessary
483 122 : IF (.NOT. ASSOCIATED(matrix_qs_ks)) THEN
484 0 : CALL dbcsr_allocate_matrix_set(matrix_qs_ks, nspin)
485 0 : DO ispin = 1, nspin
486 0 : ALLOCATE (matrix_qs_ks(ispin)%matrix)
487 : CALL dbcsr_create(matrix_qs_ks(ispin)%matrix, &
488 0 : template=matrix_qs_s(1)%matrix)
489 0 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_qs_ks(ispin)%matrix, sab_orb)
490 0 : CALL dbcsr_set(matrix_qs_ks(ispin)%matrix, 0.0_dp)
491 : END DO
492 0 : CALL set_ks_env(ks_env, matrix_ks=matrix_qs_ks)
493 : END IF
494 :
495 : ! copy to ALMO
496 250 : DO ispin = 1, nspin
497 128 : CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
498 250 : CALL dbcsr_filter(matrix_ks(ispin), eps_filter)
499 : END DO
500 :
501 122 : CALL timestop(handle)
502 :
503 122 : END SUBROUTINE init_almo_ks_matrix_via_qs
504 :
505 : ! **************************************************************************************************
506 : !> \brief Create MOs in the QS env to be able to return ALMOs to QS
507 : !> \param qs_env ...
508 : !> \param almo_scf_env ...
509 : !> \par History
510 : !> 2016.12 created [Yifei Shi]
511 : !> \author Yifei Shi
512 : ! **************************************************************************************************
513 366 : SUBROUTINE construct_qs_mos(qs_env, almo_scf_env)
514 :
515 : TYPE(qs_environment_type), POINTER :: qs_env
516 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
517 :
518 : CHARACTER(len=*), PARAMETER :: routineN = 'construct_qs_mos'
519 :
520 : INTEGER :: handle, ispin, ncol_fm, nrow_fm
521 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
522 : TYPE(cp_fm_type) :: mo_fm_copy
523 : TYPE(dft_control_type), POINTER :: dft_control
524 122 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
525 : TYPE(qs_scf_env_type), POINTER :: scf_env
526 :
527 122 : CALL timeset(routineN, handle)
528 :
529 : ! create and init scf_env (this is necessary to return MOs to qs)
530 122 : NULLIFY (mos, fm_struct_tmp, scf_env)
531 122 : ALLOCATE (scf_env)
532 122 : CALL scf_env_create(scf_env)
533 :
534 : !CALL qs_scf_env_initialize(qs_env, scf_env)
535 122 : CALL set_qs_env(qs_env, scf_env=scf_env)
536 122 : CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
537 :
538 122 : CALL dbcsr_get_info(almo_scf_env%matrix_t(1), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm)
539 :
540 : ! allocate and init mo_set
541 250 : DO ispin = 1, almo_scf_env%nspins
542 128 : CALL dbcsr_get_info(almo_scf_env%matrix_t(ispin), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm)
543 :
544 : ! Currently only fm version of mo_set is usable.
545 : ! First transform the matrix_t to fm version
546 : ! Empty the containers to prevent memory leaks
547 128 : CALL deallocate_mo_set(mos(ispin))
548 :
549 128 : IF (almo_scf_env%nspins == 1) THEN
550 : CALL allocate_mo_set(mo_set=mos(ispin), &
551 : nao=nrow_fm, &
552 : nmo=ncol_fm, &
553 : nelectron=almo_scf_env%nelectrons_total, &
554 : n_el_f=REAL(almo_scf_env%nelectrons_total, dp), &
555 : maxocc=2.0_dp, &
556 116 : flexible_electron_count=dft_control%relax_multiplicity)
557 12 : ELSE IF (almo_scf_env%nspins == 2) THEN
558 : CALL allocate_mo_set(mo_set=mos(ispin), &
559 : nao=nrow_fm, &
560 : nmo=ncol_fm, &
561 : nelectron=SUM(almo_scf_env%nocc_of_domain(:, ispin)), &
562 : n_el_f=REAL(SUM(almo_scf_env%nocc_of_domain(:, ispin)), dp), &
563 : maxocc=1.0_dp, &
564 60 : flexible_electron_count=dft_control%relax_multiplicity)
565 : END IF
566 :
567 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_fm, ncol_global=ncol_fm, &
568 : context=almo_scf_env%blacs_env, &
569 128 : para_env=almo_scf_env%para_env)
570 :
571 128 : CALL cp_fm_create(mo_fm_copy, fm_struct_tmp, name="t_orthogonal_converted_to_fm")
572 128 : CALL cp_fm_struct_release(fm_struct_tmp)
573 : !CALL copy_dbcsr_to_fm(almo_scf_env%matrix_t(ispin), mo_fm_copy)
574 :
575 128 : CALL init_mo_set(mos(ispin), fm_ref=mo_fm_copy, name='fm_mo')
576 :
577 506 : CALL cp_fm_release(mo_fm_copy)
578 :
579 : END DO
580 :
581 122 : CALL timestop(handle)
582 :
583 122 : END SUBROUTINE construct_qs_mos
584 :
585 : ! **************************************************************************************************
586 : !> \brief return density matrix to the qs_env
587 : !> \param qs_env ...
588 : !> \param matrix_p ...
589 : !> \param mat_distr_aos ...
590 : !> \par History
591 : !> 2011.05 created [Rustam Z Khaliullin]
592 : !> \author Rustam Z Khaliullin
593 : ! **************************************************************************************************
594 1778 : SUBROUTINE almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
595 : TYPE(qs_environment_type), POINTER :: qs_env
596 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
597 : INTEGER, INTENT(IN) :: mat_distr_aos
598 :
599 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_qs_env'
600 :
601 : INTEGER :: handle, ispin, nspins
602 1778 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
603 : TYPE(qs_rho_type), POINTER :: rho
604 :
605 1778 : CALL timeset(routineN, handle)
606 :
607 1778 : NULLIFY (rho, rho_ao)
608 1778 : nspins = SIZE(matrix_p)
609 1778 : CALL get_qs_env(qs_env, rho=rho)
610 1778 : CALL qs_rho_get(rho, rho_ao=rho_ao)
611 :
612 : ! set the new density matrix
613 3574 : DO ispin = 1, nspins
614 : CALL matrix_almo_to_qs(matrix_p(ispin), &
615 : rho_ao(ispin)%matrix, &
616 3574 : mat_distr_aos)
617 : END DO
618 1778 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
619 1778 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
620 :
621 1778 : CALL timestop(handle)
622 :
623 1778 : END SUBROUTINE almo_dm_to_qs_env
624 :
625 : ! **************************************************************************************************
626 : !> \brief uses the ALMO density matrix
627 : !> to compute KS matrix (inside QS environment) and the new energy
628 : !> \param qs_env ...
629 : !> \param matrix_p ...
630 : !> \param energy_total ...
631 : !> \param mat_distr_aos ...
632 : !> \param smear ...
633 : !> \param kTS_sum ...
634 : !> \par History
635 : !> 2011.05 created [Rustam Z Khaliullin]
636 : !> 2018.09 smearing support [Ruben Staub]
637 : !> \author Rustam Z Khaliullin
638 : ! **************************************************************************************************
639 1752 : SUBROUTINE almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, smear, kTS_sum)
640 : TYPE(qs_environment_type), POINTER :: qs_env
641 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
642 : REAL(KIND=dp) :: energy_total
643 : INTEGER, INTENT(IN) :: mat_distr_aos
644 : LOGICAL, INTENT(IN), OPTIONAL :: smear
645 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS_sum
646 :
647 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_qs_ks'
648 :
649 : INTEGER :: handle
650 : LOGICAL :: smearing
651 : REAL(KIND=dp) :: entropic_term
652 : TYPE(qs_energy_type), POINTER :: energy
653 :
654 1752 : CALL timeset(routineN, handle)
655 :
656 1752 : IF (PRESENT(smear)) THEN
657 1752 : smearing = smear
658 : ELSE
659 : smearing = .FALSE.
660 : END IF
661 :
662 1752 : IF (PRESENT(kTS_sum)) THEN
663 1752 : entropic_term = kTS_sum
664 : ELSE
665 : entropic_term = 0.0_dp
666 : END IF
667 :
668 1752 : NULLIFY (energy)
669 1752 : CALL get_qs_env(qs_env, energy=energy)
670 1752 : CALL almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
671 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
672 1752 : print_active=.TRUE.)
673 :
674 : !! Add electronic entropy contribution if smearing is requested
675 : !! Previous QS entropy is replaced by the sum of the entropy for each spin
676 1752 : IF (smearing) THEN
677 20 : energy%total = energy%total - energy%kTS + entropic_term
678 : END IF
679 :
680 1752 : energy_total = energy%total
681 :
682 1752 : CALL timestop(handle)
683 :
684 1752 : END SUBROUTINE almo_dm_to_qs_ks
685 :
686 : ! **************************************************************************************************
687 : !> \brief uses the ALMO density matrix
688 : !> to compute ALMO KS matrix and the new energy
689 : !> \param qs_env ...
690 : !> \param matrix_p ...
691 : !> \param matrix_ks ...
692 : !> \param energy_total ...
693 : !> \param eps_filter ...
694 : !> \param mat_distr_aos ...
695 : !> \param smear ...
696 : !> \param kTS_sum ...
697 : !> \par History
698 : !> 2011.05 created [Rustam Z Khaliullin]
699 : !> 2018.09 smearing support [Ruben Staub]
700 : !> \author Rustam Z Khaliullin
701 : ! **************************************************************************************************
702 1752 : SUBROUTINE almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, &
703 : mat_distr_aos, smear, kTS_sum)
704 :
705 : TYPE(qs_environment_type), POINTER :: qs_env
706 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p, matrix_ks
707 : REAL(KIND=dp) :: energy_total, eps_filter
708 : INTEGER, INTENT(IN) :: mat_distr_aos
709 : LOGICAL, INTENT(IN), OPTIONAL :: smear
710 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS_sum
711 :
712 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_almo_ks'
713 :
714 : INTEGER :: handle, ispin, nspins
715 : LOGICAL :: smearing
716 : REAL(KIND=dp) :: entropic_term
717 1752 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_qs_ks
718 :
719 1752 : CALL timeset(routineN, handle)
720 :
721 1752 : IF (PRESENT(smear)) THEN
722 470 : smearing = smear
723 : ELSE
724 1282 : smearing = .FALSE.
725 : END IF
726 :
727 1752 : IF (PRESENT(kTS_sum)) THEN
728 470 : entropic_term = kTS_sum
729 : ELSE
730 1282 : entropic_term = 0.0_dp
731 : END IF
732 :
733 : ! update KS matrix in the QS env
734 : CALL almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, &
735 : smear=smearing, &
736 1752 : kTS_sum=entropic_term)
737 :
738 1752 : nspins = SIZE(matrix_ks)
739 :
740 : ! get KS matrix from the QS env and convert to the ALMO format
741 1752 : CALL get_qs_env(qs_env, matrix_ks=matrix_qs_ks)
742 3522 : DO ispin = 1, nspins
743 1770 : CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
744 3522 : CALL dbcsr_filter(matrix_ks(ispin), eps_filter)
745 : END DO
746 :
747 1752 : CALL timestop(handle)
748 :
749 1752 : END SUBROUTINE almo_dm_to_almo_ks
750 :
751 : ! **************************************************************************************************
752 : !> \brief update qs_env total energy
753 : !> \param qs_env ...
754 : !> \param energy ...
755 : !> \param energy_singles_corr ...
756 : !> \par History
757 : !> 2013.03 created [Rustam Z Khaliullin]
758 : !> \author Rustam Z Khaliullin
759 : ! **************************************************************************************************
760 112 : SUBROUTINE almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr)
761 :
762 : TYPE(qs_environment_type), POINTER :: qs_env
763 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: energy, energy_singles_corr
764 :
765 : TYPE(qs_energy_type), POINTER :: qs_energy
766 :
767 112 : CALL get_qs_env(qs_env, energy=qs_energy)
768 :
769 112 : IF (PRESENT(energy_singles_corr)) THEN
770 26 : qs_energy%singles_corr = energy_singles_corr
771 : ELSE
772 86 : qs_energy%singles_corr = 0.0_dp
773 : END IF
774 :
775 112 : IF (PRESENT(energy)) THEN
776 112 : qs_energy%total = energy
777 : END IF
778 :
779 112 : qs_energy%total = qs_energy%total + qs_energy%singles_corr
780 :
781 112 : END SUBROUTINE almo_scf_update_ks_energy
782 :
783 : ! **************************************************************************************************
784 : !> \brief Creates the matrix that imposes absolute locality on MOs
785 : !> \param qs_env ...
786 : !> \param almo_scf_env ...
787 : !> \par History
788 : !> 2011.11 created [Rustam Z. Khaliullin]
789 : !> \author Rustam Z. Khaliullin
790 : ! **************************************************************************************************
791 122 : SUBROUTINE almo_scf_construct_quencher(qs_env, almo_scf_env)
792 :
793 : TYPE(qs_environment_type), POINTER :: qs_env
794 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
795 :
796 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_construct_quencher'
797 :
798 : CHARACTER :: sym
799 : INTEGER :: col, contact_atom_1, contact_atom_2, domain_col, domain_map_local_entries, &
800 : domain_row, global_entries, global_list_length, grid1, GroupID, handle, hold, iatom, &
801 : iatom2, iblock_col, iblock_row, idomain, idomain2, ientry, igrid, ineig, ineighbor, &
802 : iNode, inode2, ipair, ispin, jatom, jatom2, jdomain2, local_list_length, &
803 : max_domain_neighbors, max_neig, mynode, nblkcols_tot, nblkrows_tot, nblks, ndomains, &
804 : neig_temp, nnode2, nNodes, row, unit_nr
805 122 : INTEGER, ALLOCATABLE, DIMENSION(:) :: current_number_neighbors, domain_entries_cpu, &
806 122 : domain_map_global, domain_map_local, first_atom_of_molecule, global_list, &
807 122 : last_atom_of_molecule, list_length_cpu, list_offset_cpu, local_list, offset_for_cpu
808 122 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: domain_grid, domain_neighbor_list, &
809 122 : domain_neighbor_list_excessive
810 122 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
811 : LOGICAL :: already_listed, block_active, &
812 : delayed_increment, found, &
813 : max_neig_fails, tr
814 : REAL(KIND=dp) :: contact1_radius, contact2_radius, &
815 : distance, distance_squared, overlap, &
816 : r0, r1, s0, s1, trial_distance_squared
817 122 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: new_block
818 : REAL(KIND=dp), DIMENSION(3) :: rab
819 122 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_old_block
820 : TYPE(cell_type), POINTER :: cell
821 : TYPE(cp_logger_type), POINTER :: logger
822 : TYPE(dbcsr_distribution_type) :: dist
823 122 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
824 : TYPE(dbcsr_type) :: matrix_s_sym
825 122 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
826 : TYPE(mp_comm_type) :: group
827 : TYPE(neighbor_list_iterator_p_type), &
828 122 : DIMENSION(:), POINTER :: nl_iterator, nl_iterator2
829 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
830 122 : POINTER :: sab_almo
831 122 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
832 :
833 122 : CALL timeset(routineN, handle)
834 :
835 : ! get a useful output_unit
836 122 : logger => cp_get_default_logger()
837 122 : IF (logger%para_env%is_source()) THEN
838 61 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
839 : ELSE
840 : unit_nr = -1
841 : END IF
842 :
843 122 : ndomains = almo_scf_env%ndomains
844 :
845 : CALL get_qs_env(qs_env=qs_env, &
846 : particle_set=particle_set, &
847 : molecule_set=molecule_set, &
848 : cell=cell, &
849 : matrix_s=matrix_s, &
850 122 : sab_almo=sab_almo)
851 :
852 : ! if we are dealing with molecules get info about them
853 122 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular .OR. &
854 : almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
855 366 : ALLOCATE (first_atom_of_molecule(almo_scf_env%nmolecules))
856 244 : ALLOCATE (last_atom_of_molecule(almo_scf_env%nmolecules))
857 : CALL get_molecule_set_info(molecule_set, &
858 : mol_to_first_atom=first_atom_of_molecule, &
859 122 : mol_to_last_atom=last_atom_of_molecule)
860 : END IF
861 :
862 : ! create a symmetrized copy of the ao overlap
863 : CALL dbcsr_create(matrix_s_sym, &
864 : template=almo_scf_env%matrix_s(1), &
865 122 : matrix_type=dbcsr_type_no_symmetry)
866 : CALL dbcsr_get_info(almo_scf_env%matrix_s(1), &
867 122 : matrix_type=sym)
868 122 : IF (sym == dbcsr_type_no_symmetry) THEN
869 0 : CALL dbcsr_copy(matrix_s_sym, almo_scf_env%matrix_s(1))
870 : ELSE
871 : CALL dbcsr_desymmetrize(almo_scf_env%matrix_s(1), &
872 122 : matrix_s_sym)
873 : END IF
874 :
875 494 : ALLOCATE (almo_scf_env%quench_t(almo_scf_env%nspins))
876 494 : ALLOCATE (almo_scf_env%domain_map(almo_scf_env%nspins))
877 :
878 250 : DO ispin = 1, almo_scf_env%nspins
879 :
880 : ! create the sparsity template for the occupied orbitals
881 : CALL matrix_almo_create(matrix_new=almo_scf_env%quench_t(ispin), &
882 : matrix_qs=matrix_s(1)%matrix, &
883 : almo_scf_env=almo_scf_env, &
884 : name_new="T_QUENCHER", &
885 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_occ], &
886 : symmetry_new=dbcsr_type_no_symmetry, &
887 : spin_key=ispin, &
888 128 : init_domains=.FALSE.)
889 :
890 : ! initialize distance quencher
891 128 : CALL dbcsr_work_create(almo_scf_env%quench_t(ispin), work_mutable=.TRUE.)
892 : CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist, &
893 128 : nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
894 128 : CALL dbcsr_distribution_get(dist, numnodes=nNodes, group=GroupID, mynode=mynode)
895 128 : CALL group%set_handle(groupid)
896 :
897 : ! create global atom neighbor list from the local lists
898 : ! first, calculate number of local pairs
899 128 : local_list_length = 0
900 128 : CALL neighbor_list_iterator_create(nl_iterator, sab_almo)
901 39611 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
902 : ! nnode - total number of neighbors for iatom
903 : ! inode - current neighbor count
904 : CALL get_iterator_info(nl_iterator, &
905 39483 : iatom=iatom2, jatom=jatom2, inode=inode2, nnode=nnode2)
906 39611 : IF (inode2 == 1) THEN
907 2049 : local_list_length = local_list_length + nnode2
908 : END IF
909 : END DO
910 128 : CALL neighbor_list_iterator_release(nl_iterator)
911 :
912 : ! second, extract the local list to an array
913 383 : ALLOCATE (local_list(2*local_list_length))
914 128 : local_list(:) = 0
915 128 : local_list_length = 0
916 128 : CALL neighbor_list_iterator_create(nl_iterator2, sab_almo)
917 39611 : DO WHILE (neighbor_list_iterate(nl_iterator2) == 0)
918 : CALL get_iterator_info(nl_iterator2, &
919 39483 : iatom=iatom2, jatom=jatom2)
920 39483 : local_list(2*local_list_length + 1) = iatom2
921 39483 : local_list(2*local_list_length + 2) = jatom2
922 39483 : local_list_length = local_list_length + 1
923 : END DO ! end loop over pairs of atoms
924 128 : CALL neighbor_list_iterator_release(nl_iterator2)
925 :
926 : ! third, communicate local length to the other nodes
927 512 : ALLOCATE (list_length_cpu(nNodes), list_offset_cpu(nNodes))
928 128 : CALL group%allgather(2*local_list_length, list_length_cpu)
929 :
930 : ! fourth, create a global list
931 128 : list_offset_cpu(1) = 0
932 256 : DO iNode = 2, nNodes
933 : list_offset_cpu(iNode) = list_offset_cpu(iNode - 1) + &
934 256 : list_length_cpu(iNode - 1)
935 : END DO
936 128 : global_list_length = list_offset_cpu(nNodes) + list_length_cpu(nNodes)
937 :
938 : ! fifth, communicate all list data
939 384 : ALLOCATE (global_list(global_list_length))
940 : CALL group%allgatherv(local_list, global_list, &
941 128 : list_length_cpu, list_offset_cpu)
942 128 : DEALLOCATE (list_length_cpu, list_offset_cpu)
943 128 : DEALLOCATE (local_list)
944 :
945 : ! calculate maximum number of atoms surrounding the domain
946 384 : ALLOCATE (current_number_neighbors(almo_scf_env%ndomains))
947 128 : current_number_neighbors(:) = 0
948 128 : global_list_length = global_list_length/2
949 79094 : DO ipair = 1, global_list_length
950 78966 : iatom2 = global_list(2*(ipair - 1) + 1)
951 78966 : jatom2 = global_list(2*(ipair - 1) + 2)
952 78966 : idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
953 78966 : jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
954 : ! add to the list
955 78966 : current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
956 : ! add j,i with i,j
957 79094 : IF (idomain2 /= jdomain2) THEN
958 63144 : current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
959 : END IF
960 : END DO
961 962 : max_domain_neighbors = MAXVAL(current_number_neighbors)
962 :
963 : ! use the global atom neighbor list to create a global domain neighbor list
964 512 : ALLOCATE (domain_neighbor_list_excessive(ndomains, max_domain_neighbors))
965 962 : current_number_neighbors(:) = 1
966 962 : DO ipair = 1, ndomains
967 962 : domain_neighbor_list_excessive(ipair, 1) = ipair
968 : END DO
969 79094 : DO ipair = 1, global_list_length
970 78966 : iatom2 = global_list(2*(ipair - 1) + 1)
971 78966 : jatom2 = global_list(2*(ipair - 1) + 2)
972 78966 : idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
973 78966 : jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
974 78966 : already_listed = .FALSE.
975 325938 : DO ineighbor = 1, current_number_neighbors(idomain2)
976 325938 : IF (domain_neighbor_list_excessive(idomain2, ineighbor) == jdomain2) THEN
977 : already_listed = .TRUE.
978 : EXIT
979 : END IF
980 : END DO
981 79094 : IF (.NOT. already_listed) THEN
982 : ! add to the list
983 2722 : current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
984 2722 : domain_neighbor_list_excessive(idomain2, current_number_neighbors(idomain2)) = jdomain2
985 : ! add j,i with i,j
986 2722 : IF (idomain2 /= jdomain2) THEN
987 2722 : current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
988 2722 : domain_neighbor_list_excessive(jdomain2, current_number_neighbors(jdomain2)) = idomain2
989 : END IF
990 : END IF
991 : END DO ! end loop over pairs of atoms
992 128 : DEALLOCATE (global_list)
993 :
994 962 : max_domain_neighbors = MAXVAL(current_number_neighbors)
995 512 : ALLOCATE (domain_neighbor_list(ndomains, max_domain_neighbors))
996 128 : domain_neighbor_list(:, :) = 0
997 7248 : domain_neighbor_list(:, :) = domain_neighbor_list_excessive(:, 1:max_domain_neighbors)
998 128 : DEALLOCATE (domain_neighbor_list_excessive)
999 :
1000 384 : ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1001 384 : ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(max_domain_neighbors*ndomains, 2))
1002 12956 : almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1003 962 : almo_scf_env%domain_map(ispin)%index1(:) = 0
1004 128 : domain_map_local_entries = 0
1005 :
1006 : ! RZK-warning intermediate [0,1] quencher values are ill-defined
1007 : ! for molecules (not continuous and conceptually inadequate)
1008 :
1009 : CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), &
1010 128 : row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1011 : ! O(N) loop over domain pairs
1012 962 : DO row = 1, nblkrows_tot
1013 7240 : DO col = 1, current_number_neighbors(row)
1014 6278 : tr = .FALSE.
1015 6278 : iblock_row = row
1016 6278 : iblock_col = domain_neighbor_list(row, col)
1017 : CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin), &
1018 6278 : iblock_row, iblock_col, hold)
1019 :
1020 7112 : IF (hold == mynode) THEN
1021 :
1022 : ! Translate indices of distribution blocks to indices of domain blocks
1023 : ! Rows are AOs
1024 3139 : domain_row = almo_scf_env%domain_index_of_ao_block(iblock_row)
1025 : ! Columns are electrons (i.e. MOs)
1026 3139 : domain_col = almo_scf_env%domain_index_of_mo_block(iblock_col)
1027 :
1028 3139 : SELECT CASE (almo_scf_env%constraint_type)
1029 : CASE (almo_constraint_block_diagonal)
1030 :
1031 0 : block_active = .FALSE.
1032 : ! type of electron groups
1033 0 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1034 :
1035 : ! type of ao domains
1036 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1037 :
1038 : ! ao domains are molecular / electron groups are molecular
1039 0 : IF (domain_row == domain_col) THEN
1040 : block_active = .TRUE.
1041 : END IF
1042 :
1043 : ELSE ! ao domains are atomic
1044 :
1045 : ! ao domains are atomic / electron groups are molecular
1046 0 : CPABORT("Illegal: atomic domains and molecular groups")
1047 :
1048 : END IF
1049 :
1050 : ELSE ! electron groups are atomic
1051 :
1052 : ! type of ao domains
1053 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1054 :
1055 : ! ao domains are molecular / electron groups are atomic
1056 0 : CPABORT("Illegal: molecular domains and atomic groups")
1057 :
1058 : ELSE
1059 :
1060 : ! ao domains are atomic / electron groups are atomic
1061 0 : IF (domain_row == domain_col) THEN
1062 : block_active = .TRUE.
1063 : END IF
1064 :
1065 : END IF
1066 :
1067 : END IF ! end type of electron groups
1068 :
1069 0 : IF (block_active) THEN
1070 :
1071 0 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1072 0 : new_block(:, :) = 1.0_dp
1073 0 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1074 0 : DEALLOCATE (new_block)
1075 :
1076 0 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1077 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1078 : END IF
1079 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1080 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1081 0 : domain_map_local_entries = domain_map_local_entries + 1
1082 :
1083 : END IF
1084 :
1085 : CASE (almo_constraint_ao_overlap)
1086 :
1087 : ! type of electron groups
1088 0 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1089 :
1090 : ! type of ao domains
1091 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1092 :
1093 : ! ao domains are molecular / electron groups are molecular
1094 :
1095 : ! compute the maximum overlap between the atoms of the two molecules
1096 0 : CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1097 0 : IF (found) THEN
1098 0 : overlap = MAXVAL(ABS(p_old_block))
1099 : ELSE
1100 : overlap = 0.0_dp
1101 : END IF
1102 :
1103 : ELSE ! ao domains are atomic
1104 :
1105 : ! ao domains are atomic / electron groups are molecular
1106 : ! overlap_between_atom_and_molecule(atom=domain_row,molecule=domain_col)
1107 0 : CPABORT("atomic domains and molecular groups - NYI")
1108 :
1109 : END IF
1110 :
1111 : ELSE ! electron groups are atomic
1112 :
1113 : ! type of ao domains
1114 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1115 :
1116 : ! ao domains are molecular / electron groups are atomic
1117 : ! overlap_between_atom_and_molecule(atom=domain_col,molecule=domain_row)
1118 0 : CPABORT("molecular domains and atomic groups - NYI")
1119 :
1120 : ELSE
1121 :
1122 : ! ao domains are atomic / electron groups are atomic
1123 : ! compute max overlap between atoms: domain_row and domain_col
1124 0 : CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1125 0 : IF (found) THEN
1126 0 : overlap = MAXVAL(ABS(p_old_block))
1127 : ELSE
1128 : overlap = 0.0_dp
1129 : END IF
1130 :
1131 : END IF
1132 :
1133 : END IF ! end type of electron groups
1134 :
1135 0 : s0 = -LOG10(ABS(almo_scf_env%quencher_s0))
1136 0 : s1 = -LOG10(ABS(almo_scf_env%quencher_s1))
1137 0 : IF (overlap == 0.0_dp) THEN
1138 0 : overlap = -LOG10(ABS(almo_scf_env%eps_filter)) + 100.0_dp
1139 : ELSE
1140 0 : overlap = -LOG10(overlap)
1141 : END IF
1142 0 : IF (s0 < 0.0_dp) THEN
1143 0 : CPABORT("S0 is less than zero")
1144 : END IF
1145 0 : IF (s1 <= 0.0_dp) THEN
1146 0 : CPABORT("S1 is less than or equal to zero")
1147 : END IF
1148 0 : IF (s0 >= s1) THEN
1149 0 : CPABORT("S0 is greater than or equal to S1")
1150 : END IF
1151 :
1152 : ! Fill in non-zero blocks if AOs are close to the electron center
1153 0 : IF (overlap < s1) THEN
1154 0 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1155 0 : IF (overlap <= s0) THEN
1156 0 : new_block(:, :) = 1.0_dp
1157 : ELSE
1158 0 : new_block(:, :) = 1.0_dp/(1.0_dp + EXP(-(s0 - s1)/(s0 - overlap) - (s0 - s1)/(overlap - s1)))
1159 : END IF
1160 :
1161 0 : IF (ABS(new_block(1, 1)) > ABS(almo_scf_env%eps_filter)) THEN
1162 0 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1163 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1164 : END IF
1165 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1166 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1167 0 : domain_map_local_entries = domain_map_local_entries + 1
1168 : END IF
1169 :
1170 0 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1171 0 : DEALLOCATE (new_block)
1172 :
1173 : END IF
1174 :
1175 : CASE (almo_constraint_distance)
1176 :
1177 : ! type of electron groups
1178 3139 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1179 :
1180 : ! type of ao domains
1181 3139 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1182 :
1183 : ! ao domains are molecular / electron groups are molecular
1184 :
1185 : ! compute distance between molecules: domain_row and domain_col
1186 : ! distance between molecules is defined as the smallest
1187 : ! distance among all atom pairs
1188 3139 : IF (domain_row == domain_col) THEN
1189 417 : distance = 0.0_dp
1190 417 : contact_atom_1 = first_atom_of_molecule(domain_row)
1191 417 : contact_atom_2 = first_atom_of_molecule(domain_col)
1192 : ELSE
1193 2722 : distance_squared = 1.0E+100_dp
1194 2722 : contact_atom_1 = -1
1195 2722 : contact_atom_2 = -1
1196 9074 : DO iatom = first_atom_of_molecule(domain_row), last_atom_of_molecule(domain_row)
1197 26414 : DO jatom = first_atom_of_molecule(domain_col), last_atom_of_molecule(domain_col)
1198 17340 : rab(:) = pbc(particle_set(iatom)%r(:), particle_set(jatom)%r(:), cell)
1199 17340 : trial_distance_squared = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1200 23692 : IF (trial_distance_squared < distance_squared) THEN
1201 6391 : distance_squared = trial_distance_squared
1202 6391 : contact_atom_1 = iatom
1203 6391 : contact_atom_2 = jatom
1204 : END IF
1205 : END DO ! jatom
1206 : END DO ! iatom
1207 2722 : CPASSERT(contact_atom_1 > 0)
1208 2722 : distance = SQRT(distance_squared)
1209 : END IF
1210 :
1211 : ELSE ! ao domains are atomic
1212 :
1213 : ! ao domains are atomic / electron groups are molecular
1214 : !distance_between_atom_and_molecule(atom=domain_row,molecule=domain_col)
1215 0 : CPABORT("atomic domains and molecular groups - NYI")
1216 :
1217 : END IF
1218 :
1219 : ELSE ! electron groups are atomic
1220 :
1221 : ! type of ao domains
1222 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1223 :
1224 : ! ao domains are molecular / electron groups are atomic
1225 : !distance_between_atom_and_molecule(atom=domain_col,molecule=domain_row)
1226 0 : CPABORT("molecular domains and atomic groups - NYI")
1227 :
1228 : ELSE
1229 :
1230 : ! ao domains are atomic / electron groups are atomic
1231 : ! compute distance between atoms: domain_row and domain_col
1232 0 : rab(:) = pbc(particle_set(domain_row)%r(:), particle_set(domain_col)%r(:), cell)
1233 0 : distance = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1234 0 : contact_atom_1 = domain_row
1235 0 : contact_atom_2 = domain_col
1236 :
1237 : END IF
1238 :
1239 : END IF ! end type of electron groups
1240 :
1241 : ! get atomic radii to compute distance cutoff threshold
1242 3139 : IF (almo_scf_env%quencher_radius_type == do_bondparm_covalent) THEN
1243 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1244 0 : rcov=contact1_radius)
1245 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1246 0 : rcov=contact2_radius)
1247 3139 : ELSE IF (almo_scf_env%quencher_radius_type == do_bondparm_vdw) THEN
1248 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1249 3139 : rvdw=contact1_radius)
1250 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1251 3139 : rvdw=contact2_radius)
1252 : ELSE
1253 0 : CPABORT("Illegal quencher_radius_type")
1254 : END IF
1255 3139 : contact1_radius = cp_unit_to_cp2k(contact1_radius, "angstrom")
1256 3139 : contact2_radius = cp_unit_to_cp2k(contact2_radius, "angstrom")
1257 :
1258 : !RZK-warning the procedure is faulty for molecules:
1259 : ! the closest contacts should be found using
1260 : ! the element specific radii
1261 :
1262 : ! compute inner and outer cutoff radii
1263 3139 : r0 = almo_scf_env%quencher_r0_factor*(contact1_radius + contact2_radius)
1264 : !+almo_scf_env%quencher_r0_shift
1265 3139 : r1 = almo_scf_env%quencher_r1_factor*(contact1_radius + contact2_radius)
1266 : !+almo_scf_env%quencher_r1_shift
1267 :
1268 3139 : IF (r0 < 0.0_dp) THEN
1269 0 : CPABORT("R0 is less than zero")
1270 : END IF
1271 3139 : IF (r1 <= 0.0_dp) THEN
1272 0 : CPABORT("R1 is less than or equal to zero")
1273 : END IF
1274 3139 : IF (r0 > r1) THEN
1275 0 : CPABORT("R0 is greater than or equal to R1")
1276 : END IF
1277 :
1278 : ! Fill in non-zero blocks if AOs are close to the electron center
1279 3139 : IF (distance < r1) THEN
1280 8740 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1281 2185 : IF (distance <= r0) THEN
1282 101919 : new_block(:, :) = 1.0_dp
1283 : ELSE
1284 : ! remove the intermediate values from the quencher temporarily
1285 0 : CPABORT("distance > r0 not yet validated") ! Unexplained in https://github.com/cp2k/cp2k/pull/5345
1286 0 : new_block(:, :) = 1.0_dp/(1.0_dp + EXP((r1 - r0)/(r0 - distance) + (r1 - r0)/(r1 - distance)))
1287 : END IF
1288 :
1289 2185 : IF (ABS(new_block(1, 1)) > ABS(almo_scf_env%eps_filter)) THEN
1290 2185 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1291 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1292 : END IF
1293 2185 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1294 2185 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1295 2185 : domain_map_local_entries = domain_map_local_entries + 1
1296 : END IF
1297 :
1298 2185 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1299 2185 : DEALLOCATE (new_block)
1300 : END IF
1301 :
1302 : CASE DEFAULT
1303 3139 : CPABORT("Illegal constraint type")
1304 : END SELECT
1305 :
1306 : END IF ! mynode
1307 :
1308 : END DO
1309 : END DO ! end O(N) loop over pairs
1310 :
1311 128 : DEALLOCATE (domain_neighbor_list)
1312 128 : DEALLOCATE (current_number_neighbors)
1313 :
1314 128 : CALL dbcsr_finalize(almo_scf_env%quench_t(ispin))
1315 :
1316 : CALL dbcsr_filter(almo_scf_env%quench_t(ispin), &
1317 128 : almo_scf_env%eps_filter)
1318 :
1319 : ! check that both domain_map and quench_t have the same number of entries
1320 128 : nblks = dbcsr_get_num_blocks(almo_scf_env%quench_t(ispin))
1321 128 : IF (nblks /= domain_map_local_entries) THEN
1322 0 : CPABORT("number of blocks is wrong")
1323 : END IF
1324 :
1325 : ! first, communicate map sizes on the other nodes
1326 384 : ALLOCATE (domain_entries_cpu(nNodes), offset_for_cpu(nNodes))
1327 128 : CALL group%allgather(2*domain_map_local_entries, domain_entries_cpu)
1328 :
1329 : ! second, create
1330 128 : offset_for_cpu(1) = 0
1331 256 : DO iNode = 2, nNodes
1332 : offset_for_cpu(iNode) = offset_for_cpu(iNode - 1) + &
1333 256 : domain_entries_cpu(iNode - 1)
1334 : END DO
1335 128 : global_entries = offset_for_cpu(nNodes) + domain_entries_cpu(nNodes)
1336 :
1337 : ! communicate all entries
1338 384 : ALLOCATE (domain_map_global(global_entries))
1339 508 : ALLOCATE (domain_map_local(2*domain_map_local_entries))
1340 2313 : DO ientry = 1, domain_map_local_entries
1341 2185 : domain_map_local(2*(ientry - 1) + 1) = almo_scf_env%domain_map(ispin)%pairs(ientry, 1)
1342 2313 : domain_map_local(2*ientry) = almo_scf_env%domain_map(ispin)%pairs(ientry, 2)
1343 : END DO
1344 : CALL group%allgatherv(domain_map_local, domain_map_global, &
1345 128 : domain_entries_cpu, offset_for_cpu)
1346 128 : DEALLOCATE (domain_entries_cpu, offset_for_cpu)
1347 128 : DEALLOCATE (domain_map_local)
1348 :
1349 128 : DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
1350 128 : DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
1351 256 : ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1352 512 : ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(global_entries/2, 2))
1353 9124 : almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1354 962 : almo_scf_env%domain_map(ispin)%index1(:) = 0
1355 :
1356 : ! unpack the received data into a local variable
1357 : ! since we do not know the maximum global number of neighbors
1358 : ! try one. if fails increase the maximum number and try again
1359 : ! until it succeeds
1360 : max_neig = max_domain_neighbors
1361 : max_neig_fails = .TRUE.
1362 256 : max_neig_loop: DO WHILE (max_neig_fails)
1363 512 : ALLOCATE (domain_grid(almo_scf_env%ndomains, 0:max_neig))
1364 128 : domain_grid(:, :) = 0
1365 : ! init the number of collected neighbors
1366 962 : domain_grid(:, 0) = 1
1367 : ! loop over the records
1368 128 : global_entries = global_entries/2
1369 4498 : DO ientry = 1, global_entries
1370 : ! get the center
1371 4370 : grid1 = domain_map_global(2*ientry)
1372 : ! get the neighbor
1373 4370 : ineig = domain_map_global(2*(ientry - 1) + 1)
1374 : ! check boundaries
1375 4370 : IF (domain_grid(grid1, 0) > max_neig) THEN
1376 : ! this neighbor will overstep the boundaries
1377 : ! stop the trial and increase the max number of neighbors
1378 0 : DEALLOCATE (domain_grid)
1379 0 : max_neig = max_neig*2
1380 0 : CYCLE max_neig_loop
1381 : END IF
1382 : ! for the current center loop over the collected neighbors
1383 : ! to insert the current record in a numerical order
1384 : delayed_increment = .FALSE.
1385 19016 : DO igrid = 1, domain_grid(grid1, 0)
1386 : ! compare the current neighbor with that already in the 'book'
1387 19016 : IF (ineig < domain_grid(grid1, igrid)) THEN
1388 : ! if this one is smaller then insert it here and pick up the one
1389 : ! from the book to continue inserting
1390 3180 : neig_temp = ineig
1391 3180 : ineig = domain_grid(grid1, igrid)
1392 3180 : domain_grid(grid1, igrid) = neig_temp
1393 : ELSE
1394 11466 : IF (domain_grid(grid1, igrid) == 0) THEN
1395 : ! got the empty slot now - insert the record
1396 4370 : domain_grid(grid1, igrid) = ineig
1397 : ! increase the record counter but do it outside the loop
1398 4370 : delayed_increment = .TRUE.
1399 : END IF
1400 : END IF
1401 : END DO
1402 4498 : IF (delayed_increment) THEN
1403 4370 : domain_grid(grid1, 0) = domain_grid(grid1, 0) + 1
1404 : ELSE
1405 : ! should not be here - all records must be inserted
1406 0 : CPABORT("all records must be inserted")
1407 : END IF
1408 : END DO
1409 128 : max_neig_fails = .FALSE.
1410 : END DO max_neig_loop
1411 128 : DEALLOCATE (domain_map_global)
1412 :
1413 128 : ientry = 1
1414 962 : DO idomain = 1, almo_scf_env%ndomains
1415 5204 : DO ineig = 1, domain_grid(idomain, 0) - 1
1416 4370 : almo_scf_env%domain_map(ispin)%pairs(ientry, 1) = domain_grid(idomain, ineig)
1417 4370 : almo_scf_env%domain_map(ispin)%pairs(ientry, 2) = idomain
1418 5204 : ientry = ientry + 1
1419 : END DO
1420 962 : almo_scf_env%domain_map(ispin)%index1(idomain) = ientry
1421 : END DO
1422 378 : DEALLOCATE (domain_grid)
1423 :
1424 : END DO ! ispin
1425 122 : IF (almo_scf_env%nspins == 2) THEN
1426 : CALL dbcsr_copy(almo_scf_env%quench_t(2), &
1427 6 : almo_scf_env%quench_t(1))
1428 : almo_scf_env%domain_map(2)%pairs(:, :) = &
1429 50 : almo_scf_env%domain_map(1)%pairs(:, :)
1430 : almo_scf_env%domain_map(2)%index1(:) = &
1431 18 : almo_scf_env%domain_map(1)%index1(:)
1432 : END IF
1433 :
1434 122 : CALL dbcsr_release(matrix_s_sym)
1435 :
1436 122 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular .OR. &
1437 : almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1438 122 : DEALLOCATE (first_atom_of_molecule)
1439 122 : DEALLOCATE (last_atom_of_molecule)
1440 : END IF
1441 :
1442 122 : CALL timestop(handle)
1443 :
1444 244 : END SUBROUTINE almo_scf_construct_quencher
1445 :
1446 : ! *****************************************************************************
1447 : !> \brief Compute matrix W (energy-weighted density matrix) that is needed
1448 : !> for the evaluation of forces
1449 : !> \param matrix_w ...
1450 : !> \param almo_scf_env ...
1451 : !> \par History
1452 : !> 2015.03 created [Rustam Z. Khaliullin]
1453 : !> \author Rustam Z. Khaliullin
1454 : ! **************************************************************************************************
1455 66 : SUBROUTINE calculate_w_matrix_almo(matrix_w, almo_scf_env)
1456 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
1457 : TYPE(almo_scf_env_type) :: almo_scf_env
1458 :
1459 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_almo'
1460 :
1461 : INTEGER :: handle, ispin
1462 : REAL(KIND=dp) :: scaling
1463 : TYPE(dbcsr_type) :: tmp_nn1, tmp_no1, tmp_oo1, tmp_oo2
1464 :
1465 66 : CALL timeset(routineN, handle)
1466 :
1467 66 : IF (almo_scf_env%nspins == 1) THEN
1468 66 : scaling = 2.0_dp
1469 : ELSE
1470 0 : scaling = 1.0_dp
1471 : END IF
1472 :
1473 132 : DO ispin = 1, almo_scf_env%nspins
1474 :
1475 : CALL dbcsr_create(tmp_nn1, template=almo_scf_env%matrix_s(1), &
1476 66 : matrix_type=dbcsr_type_no_symmetry)
1477 : CALL dbcsr_create(tmp_no1, template=almo_scf_env%matrix_t(ispin), &
1478 66 : matrix_type=dbcsr_type_no_symmetry)
1479 : CALL dbcsr_create(tmp_oo1, template=almo_scf_env%matrix_sigma_inv(ispin), &
1480 66 : matrix_type=dbcsr_type_no_symmetry)
1481 : CALL dbcsr_create(tmp_oo2, template=almo_scf_env%matrix_sigma_inv(ispin), &
1482 66 : matrix_type=dbcsr_type_no_symmetry)
1483 :
1484 66 : CALL dbcsr_copy(tmp_nn1, almo_scf_env%matrix_ks(ispin))
1485 : ! 1. TMP_NO1=F.T
1486 : CALL dbcsr_multiply("N", "N", scaling, tmp_nn1, almo_scf_env%matrix_t(ispin), &
1487 66 : 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1488 : ! 2. TMP_OO1=T^(tr).TMP_NO1=T^(tr).(FT)
1489 : CALL dbcsr_multiply("T", "N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_no1, &
1490 66 : 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1491 : ! 3. TMP_OO2=TMP_OO1.siginv=(T^(tr)FT).siginv
1492 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_oo1, almo_scf_env%matrix_sigma_inv(ispin), &
1493 66 : 0.0_dp, tmp_oo2, filter_eps=almo_scf_env%eps_filter)
1494 : ! 4. TMP_OO1=siginv.TMP_OO2=siginv.(T^(tr)FTsiginv)
1495 : CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_sigma_inv(ispin), tmp_oo2, &
1496 66 : 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1497 : ! 5. TMP_NO1=T.TMP_OO1.=T.(siginvT^(tr)FTsiginv)
1498 : CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_oo1, &
1499 66 : 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1500 : ! 6. TMP_NN1=TMP_NO1.T^(tr)=(TsiginvT^(tr)FTsiginv).T^(tr)=RFR
1501 : CALL dbcsr_multiply("N", "T", 1.0_dp, tmp_no1, almo_scf_env%matrix_t(ispin), &
1502 66 : 0.0_dp, tmp_nn1, filter_eps=almo_scf_env%eps_filter)
1503 66 : CALL matrix_almo_to_qs(tmp_nn1, matrix_w(ispin)%matrix, almo_scf_env%mat_distr_aos)
1504 :
1505 66 : CALL dbcsr_release(tmp_nn1)
1506 66 : CALL dbcsr_release(tmp_no1)
1507 66 : CALL dbcsr_release(tmp_oo1)
1508 132 : CALL dbcsr_release(tmp_oo2)
1509 :
1510 : END DO
1511 :
1512 66 : CALL timestop(handle)
1513 :
1514 66 : END SUBROUTINE calculate_w_matrix_almo
1515 :
1516 : END MODULE almo_scf_qs
1517 :
|
