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 : ! startQQQ - this part of the code scales quadratically
292 : ! therefore it is replaced with a less general but linear scaling algorithm below
293 : ! the quadratic algorithm is kept to be re-written later
294 : !QQQCALL dbcsr_get_info(matrix_new, nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
295 : !QQQDO row = 1, nblkrows_tot
296 : !QQQ DO col = 1, nblkcols_tot
297 : !QQQ tr = .FALSE.
298 : !QQQ iblock_row = row
299 : !QQQ iblock_col = col
300 : !QQQ CALL dbcsr_get_stored_coordinates(matrix_new, iblock_row, iblock_col, tr, hold)
301 :
302 : !QQQ IF(hold==mynode) THEN
303 : !QQQ
304 : !QQQ ! RZK-warning replace with a function which says if this
305 : !QQQ ! distribution block is active or not
306 : !QQQ ! Translate indeces of distribution blocks to domain blocks
307 : !QQQ if (size_keys(1)==almo_mat_dim_aobasis) then
308 : !QQQ domain_row=almo_scf_env%domain_index_of_ao_block(iblock_row)
309 : !QQQ else if (size_keys(2)==almo_mat_dim_occ .OR. &
310 : !QQQ size_keys(2)==almo_mat_dim_virt .OR. &
311 : !QQQ size_keys(2)==almo_mat_dim_virt_disc .OR. &
312 : !QQQ size_keys(2)==almo_mat_dim_virt_full) then
313 : !QQQ domain_row=almo_scf_env%domain_index_of_mo_block(iblock_row)
314 : !QQQ else
315 : !QQQ CPErrorMessage(cp_failure_level,routineP,"Illegal dimension")
316 : !QQQ CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
317 : !QQQ endif
318 :
319 : !QQQ if (size_keys(2)==almo_mat_dim_aobasis) then
320 : !QQQ domain_col=almo_scf_env%domain_index_of_ao_block(iblock_col)
321 : !QQQ else if (size_keys(2)==almo_mat_dim_occ .OR. &
322 : !QQQ size_keys(2)==almo_mat_dim_virt .OR. &
323 : !QQQ size_keys(2)==almo_mat_dim_virt_disc .OR. &
324 : !QQQ size_keys(2)==almo_mat_dim_virt_full) then
325 : !QQQ domain_col=almo_scf_env%domain_index_of_mo_block(iblock_col)
326 : !QQQ else
327 : !QQQ CPErrorMessage(cp_failure_level,routineP,"Illegal dimension")
328 : !QQQ CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
329 : !QQQ endif
330 :
331 : !QQQ ! Finds if we need this block
332 : !QQQ ! only the block-diagonal constraint is implemented here
333 : !QQQ active=.false.
334 : !QQQ if (domain_row==domain_col) active=.true.
335 :
336 : !QQQ IF (active) THEN
337 : !QQQ ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
338 : !QQQ new_block(:, :) = 1.0_dp
339 : !QQQ CALL dbcsr_put_block(matrix_new, iblock_row, iblock_col, new_block)
340 : !QQQ DEALLOCATE (new_block)
341 : !QQQ ENDIF
342 :
343 : !QQQ ENDIF ! mynode
344 : !QQQ ENDDO
345 : !QQQENDDO
346 : !QQQtake care of zero-electron fragments
347 : ! endQQQ - end of the quadratic part
348 : ! start linear-scaling replacement:
349 : ! works only for molecular blocks AND molecular distributions
350 10870 : DO row = 1, nblkrows_tot
351 9444 : tr = .FALSE.
352 9444 : iblock_row = row
353 9444 : iblock_col = row
354 9444 : CALL dbcsr_get_stored_coordinates(matrix_new, iblock_row, iblock_col, hold)
355 :
356 10870 : IF (hold == mynode) THEN
357 :
358 14166 : active = .TRUE.
359 :
360 : one_dim_is_mo = .FALSE.
361 14166 : DO dimen = 1, 2 ! 1 - row, 2 - column dimension
362 14166 : IF (size_keys(dimen) == almo_mat_dim_occ) one_dim_is_mo = .TRUE.
363 : END DO
364 4722 : IF (one_dim_is_mo) THEN
365 1668 : IF (almo_scf_env%nocc_of_domain(row, spin_key) == 0) active = .FALSE.
366 : END IF
367 :
368 4722 : one_dim_is_mo = .FALSE.
369 14166 : DO dimen = 1, 2
370 14166 : IF (size_keys(dimen) == almo_mat_dim_virt) one_dim_is_mo = .TRUE.
371 : END DO
372 4722 : IF (one_dim_is_mo) THEN
373 417 : IF (almo_scf_env%nvirt_of_domain(row, spin_key) == 0) active = .FALSE.
374 : END IF
375 :
376 4722 : one_dim_is_mo = .FALSE.
377 14166 : DO dimen = 1, 2
378 14166 : IF (size_keys(dimen) == almo_mat_dim_virt_disc) one_dim_is_mo = .TRUE.
379 : END DO
380 4722 : IF (one_dim_is_mo) THEN
381 0 : IF (almo_scf_env%nvirt_disc_of_domain(row, spin_key) == 0) active = .FALSE.
382 : END IF
383 :
384 4722 : one_dim_is_mo = .FALSE.
385 14166 : DO dimen = 1, 2
386 14166 : IF (size_keys(dimen) == almo_mat_dim_virt_full) one_dim_is_mo = .TRUE.
387 : END DO
388 4722 : IF (one_dim_is_mo) THEN
389 417 : IF (almo_scf_env%nvirt_full_of_domain(row, spin_key) == 0) active = .FALSE.
390 : END IF
391 :
392 4688 : IF (active) THEN
393 17592 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
394 840230 : new_block(:, :) = 1.0_dp
395 4398 : CALL dbcsr_put_block(matrix_new, iblock_row, iblock_col, new_block)
396 4398 : DEALLOCATE (new_block)
397 : END IF
398 :
399 : END IF ! mynode
400 : END DO
401 : ! end lnear-scaling replacement
402 :
403 : END IF ! init_domains
404 :
405 3596 : CALL dbcsr_finalize(matrix_new)
406 :
407 3596 : CALL timestop(handle)
408 :
409 7192 : END SUBROUTINE matrix_almo_create
410 :
411 : ! **************************************************************************************************
412 : !> \brief convert between two types of matrices: QS style to ALMO style
413 : !> \param matrix_qs ...
414 : !> \param matrix_almo ...
415 : !> \param mat_distr_aos ...
416 : !> \par History
417 : !> 2011.06 created [Rustam Z Khaliullin]
418 : !> \author Rustam Z Khaliullin
419 : ! **************************************************************************************************
420 2080 : SUBROUTINE matrix_qs_to_almo(matrix_qs, matrix_almo, mat_distr_aos)
421 :
422 : TYPE(dbcsr_type) :: matrix_qs, matrix_almo
423 : INTEGER :: mat_distr_aos
424 :
425 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_qs_to_almo'
426 :
427 : INTEGER :: handle
428 : TYPE(dbcsr_type) :: matrix_qs_nosym
429 :
430 2080 : CALL timeset(routineN, handle)
431 : !RZK-warning if it's not a N(AO)xN(AO) matrix then stop
432 :
433 2080 : SELECT CASE (mat_distr_aos)
434 : CASE (almo_mat_distr_atomic)
435 : ! automatic data_type conversion
436 0 : CALL dbcsr_copy(matrix_almo, matrix_qs)
437 : CASE (almo_mat_distr_molecular)
438 : ! desymmetrize the qs matrix
439 2080 : CALL dbcsr_create(matrix_qs_nosym, template=matrix_qs, matrix_type=dbcsr_type_no_symmetry)
440 2080 : CALL dbcsr_desymmetrize(matrix_qs, matrix_qs_nosym)
441 :
442 : ! perform the magic complete_redistribute
443 : ! before calling complete_redistribute set all blocks to zero
444 : ! otherwise the non-zero elements of the redistributed matrix,
445 : ! which are in zero-blocks of the original matrix, will remain
446 : ! in the final redistributed matrix. this is a bug in
447 : ! complete_redistribute. RZK-warning it should be later corrected by calling
448 : ! dbcsr_set to 0.0 from within complete_redistribute
449 2080 : CALL dbcsr_set(matrix_almo, 0.0_dp)
450 2080 : CALL dbcsr_complete_redistribute(matrix_qs_nosym, matrix_almo)
451 2080 : CALL dbcsr_release(matrix_qs_nosym)
452 :
453 : CASE DEFAULT
454 2080 : CPABORT("")
455 : END SELECT
456 :
457 2080 : CALL timestop(handle)
458 :
459 2080 : END SUBROUTINE matrix_qs_to_almo
460 :
461 : ! **************************************************************************************************
462 : !> \brief convert between two types of matrices: ALMO style to QS style
463 : !> \param matrix_almo ...
464 : !> \param matrix_qs ...
465 : !> \param mat_distr_aos ...
466 : !> \par History
467 : !> 2011.06 created [Rustam Z Khaliullin]
468 : !> \author Rustam Z Khaliullin
469 : ! **************************************************************************************************
470 1862 : SUBROUTINE matrix_almo_to_qs(matrix_almo, matrix_qs, mat_distr_aos)
471 : TYPE(dbcsr_type) :: matrix_almo, matrix_qs
472 : INTEGER, INTENT(IN) :: mat_distr_aos
473 :
474 : CHARACTER(len=*), PARAMETER :: routineN = 'matrix_almo_to_qs'
475 :
476 : INTEGER :: handle
477 : TYPE(dbcsr_type) :: matrix_almo_redist
478 :
479 1862 : CALL timeset(routineN, handle)
480 : ! RZK-warning if it's not a N(AO)xN(AO) matrix then stop
481 :
482 1862 : SELECT CASE (mat_distr_aos)
483 : CASE (almo_mat_distr_atomic)
484 0 : CALL dbcsr_copy(matrix_qs, matrix_almo, keep_sparsity=.TRUE.)
485 : CASE (almo_mat_distr_molecular)
486 1862 : CALL dbcsr_create(matrix_almo_redist, template=matrix_qs)
487 1862 : CALL dbcsr_complete_redistribute(matrix_almo, matrix_almo_redist)
488 1862 : CALL dbcsr_set(matrix_qs, 0.0_dp)
489 1862 : CALL dbcsr_copy(matrix_qs, matrix_almo_redist, keep_sparsity=.TRUE.)
490 1862 : CALL dbcsr_release(matrix_almo_redist)
491 : CASE DEFAULT
492 1862 : CPABORT("")
493 : END SELECT
494 :
495 1862 : CALL timestop(handle)
496 :
497 1862 : END SUBROUTINE matrix_almo_to_qs
498 :
499 : ! **************************************************************************************************
500 : !> \brief Initialization of the QS and ALMO KS matrix
501 : !> \param qs_env ...
502 : !> \param matrix_ks ...
503 : !> \param mat_distr_aos ...
504 : !> \param eps_filter ...
505 : !> \par History
506 : !> 2011.05 created [Rustam Z Khaliullin]
507 : !> \author Rustam Z Khaliullin
508 : ! **************************************************************************************************
509 122 : SUBROUTINE init_almo_ks_matrix_via_qs(qs_env, matrix_ks, mat_distr_aos, eps_filter)
510 :
511 : TYPE(qs_environment_type), POINTER :: qs_env
512 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_ks
513 : INTEGER :: mat_distr_aos
514 : REAL(KIND=dp) :: eps_filter
515 :
516 : CHARACTER(len=*), PARAMETER :: routineN = 'init_almo_ks_matrix_via_qs'
517 :
518 : INTEGER :: handle, ispin, nspin
519 122 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_qs_ks, matrix_qs_s
520 : TYPE(dft_control_type), POINTER :: dft_control
521 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
522 122 : POINTER :: sab_orb
523 : TYPE(qs_ks_env_type), POINTER :: ks_env
524 :
525 122 : CALL timeset(routineN, handle)
526 :
527 122 : NULLIFY (sab_orb)
528 :
529 : ! get basic quantities from the qs_env
530 : CALL get_qs_env(qs_env, &
531 : dft_control=dft_control, &
532 : matrix_s=matrix_qs_s, &
533 : matrix_ks=matrix_qs_ks, &
534 : ks_env=ks_env, &
535 122 : sab_orb=sab_orb)
536 :
537 122 : nspin = dft_control%nspins
538 :
539 : ! create matrix_ks in the QS env if necessary
540 122 : IF (.NOT. ASSOCIATED(matrix_qs_ks)) THEN
541 0 : CALL dbcsr_allocate_matrix_set(matrix_qs_ks, nspin)
542 0 : DO ispin = 1, nspin
543 0 : ALLOCATE (matrix_qs_ks(ispin)%matrix)
544 : CALL dbcsr_create(matrix_qs_ks(ispin)%matrix, &
545 0 : template=matrix_qs_s(1)%matrix)
546 0 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_qs_ks(ispin)%matrix, sab_orb)
547 0 : CALL dbcsr_set(matrix_qs_ks(ispin)%matrix, 0.0_dp)
548 : END DO
549 0 : CALL set_ks_env(ks_env, matrix_ks=matrix_qs_ks)
550 : END IF
551 :
552 : ! copy to ALMO
553 250 : DO ispin = 1, nspin
554 128 : CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
555 250 : CALL dbcsr_filter(matrix_ks(ispin), eps_filter)
556 : END DO
557 :
558 122 : CALL timestop(handle)
559 :
560 122 : END SUBROUTINE init_almo_ks_matrix_via_qs
561 :
562 : ! **************************************************************************************************
563 : !> \brief Create MOs in the QS env to be able to return ALMOs to QS
564 : !> \param qs_env ...
565 : !> \param almo_scf_env ...
566 : !> \par History
567 : !> 2016.12 created [Yifei Shi]
568 : !> \author Yifei Shi
569 : ! **************************************************************************************************
570 366 : SUBROUTINE construct_qs_mos(qs_env, almo_scf_env)
571 :
572 : TYPE(qs_environment_type), POINTER :: qs_env
573 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
574 :
575 : CHARACTER(len=*), PARAMETER :: routineN = 'construct_qs_mos'
576 :
577 : INTEGER :: handle, ispin, ncol_fm, nrow_fm
578 : TYPE(cp_fm_struct_type), POINTER :: fm_struct_tmp
579 : TYPE(cp_fm_type) :: mo_fm_copy
580 : TYPE(dft_control_type), POINTER :: dft_control
581 122 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
582 : TYPE(qs_scf_env_type), POINTER :: scf_env
583 :
584 122 : CALL timeset(routineN, handle)
585 :
586 : ! create and init scf_env (this is necessary to return MOs to qs)
587 122 : NULLIFY (mos, fm_struct_tmp, scf_env)
588 122 : ALLOCATE (scf_env)
589 122 : CALL scf_env_create(scf_env)
590 :
591 : !CALL qs_scf_env_initialize(qs_env, scf_env)
592 122 : CALL set_qs_env(qs_env, scf_env=scf_env)
593 122 : CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
594 :
595 122 : CALL dbcsr_get_info(almo_scf_env%matrix_t(1), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm)
596 :
597 : ! allocate and init mo_set
598 250 : DO ispin = 1, almo_scf_env%nspins
599 128 : CALL dbcsr_get_info(almo_scf_env%matrix_t(ispin), nfullrows_total=nrow_fm, nfullcols_total=ncol_fm)
600 :
601 : ! Currently only fm version of mo_set is usable.
602 : ! First transform the matrix_t to fm version
603 : ! Empty the containers to prevent memory leaks
604 128 : CALL deallocate_mo_set(mos(ispin))
605 :
606 128 : IF (almo_scf_env%nspins == 1) THEN
607 : CALL allocate_mo_set(mo_set=mos(ispin), &
608 : nao=nrow_fm, &
609 : nmo=ncol_fm, &
610 : nelectron=almo_scf_env%nelectrons_total, &
611 : n_el_f=REAL(almo_scf_env%nelectrons_total, dp), &
612 : maxocc=2.0_dp, &
613 116 : flexible_electron_count=dft_control%relax_multiplicity)
614 12 : ELSEIF (almo_scf_env%nspins == 2) THEN
615 : CALL allocate_mo_set(mo_set=mos(ispin), &
616 : nao=nrow_fm, &
617 : nmo=ncol_fm, &
618 : nelectron=SUM(almo_scf_env%nocc_of_domain(:, ispin)), &
619 : n_el_f=REAL(SUM(almo_scf_env%nocc_of_domain(:, ispin)), dp), &
620 : maxocc=1.0_dp, &
621 60 : flexible_electron_count=dft_control%relax_multiplicity)
622 : END IF
623 :
624 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_fm, ncol_global=ncol_fm, &
625 : context=almo_scf_env%blacs_env, &
626 128 : para_env=almo_scf_env%para_env)
627 :
628 128 : CALL cp_fm_create(mo_fm_copy, fm_struct_tmp, name="t_orthogonal_converted_to_fm")
629 128 : CALL cp_fm_struct_release(fm_struct_tmp)
630 : !CALL copy_dbcsr_to_fm(almo_scf_env%matrix_t(ispin), mo_fm_copy)
631 :
632 128 : CALL init_mo_set(mos(ispin), fm_ref=mo_fm_copy, name='fm_mo')
633 :
634 506 : CALL cp_fm_release(mo_fm_copy)
635 :
636 : END DO
637 :
638 122 : CALL timestop(handle)
639 :
640 122 : END SUBROUTINE construct_qs_mos
641 :
642 : ! **************************************************************************************************
643 : !> \brief return density matrix to the qs_env
644 : !> \param qs_env ...
645 : !> \param matrix_p ...
646 : !> \param mat_distr_aos ...
647 : !> \par History
648 : !> 2011.05 created [Rustam Z Khaliullin]
649 : !> \author Rustam Z Khaliullin
650 : ! **************************************************************************************************
651 1778 : SUBROUTINE almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
652 : TYPE(qs_environment_type), POINTER :: qs_env
653 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
654 : INTEGER, INTENT(IN) :: mat_distr_aos
655 :
656 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_qs_env'
657 :
658 : INTEGER :: handle, ispin, nspins
659 1778 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
660 : TYPE(qs_rho_type), POINTER :: rho
661 :
662 1778 : CALL timeset(routineN, handle)
663 :
664 1778 : NULLIFY (rho, rho_ao)
665 1778 : nspins = SIZE(matrix_p)
666 1778 : CALL get_qs_env(qs_env, rho=rho)
667 1778 : CALL qs_rho_get(rho, rho_ao=rho_ao)
668 :
669 : ! set the new density matrix
670 3574 : DO ispin = 1, nspins
671 : CALL matrix_almo_to_qs(matrix_p(ispin), &
672 : rho_ao(ispin)%matrix, &
673 3574 : mat_distr_aos)
674 : END DO
675 1778 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
676 1778 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
677 :
678 1778 : CALL timestop(handle)
679 :
680 1778 : END SUBROUTINE almo_dm_to_qs_env
681 :
682 : ! **************************************************************************************************
683 : !> \brief uses the ALMO density matrix
684 : !> to compute KS matrix (inside QS environment) and the new energy
685 : !> \param qs_env ...
686 : !> \param matrix_p ...
687 : !> \param energy_total ...
688 : !> \param mat_distr_aos ...
689 : !> \param smear ...
690 : !> \param kTS_sum ...
691 : !> \par History
692 : !> 2011.05 created [Rustam Z Khaliullin]
693 : !> 2018.09 smearing support [Ruben Staub]
694 : !> \author Rustam Z Khaliullin
695 : ! **************************************************************************************************
696 1752 : SUBROUTINE almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, smear, kTS_sum)
697 : TYPE(qs_environment_type), POINTER :: qs_env
698 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
699 : REAL(KIND=dp) :: energy_total
700 : INTEGER, INTENT(IN) :: mat_distr_aos
701 : LOGICAL, INTENT(IN), OPTIONAL :: smear
702 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS_sum
703 :
704 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_qs_ks'
705 :
706 : INTEGER :: handle
707 : LOGICAL :: smearing
708 : REAL(KIND=dp) :: entropic_term
709 : TYPE(qs_energy_type), POINTER :: energy
710 :
711 1752 : CALL timeset(routineN, handle)
712 :
713 1752 : IF (PRESENT(smear)) THEN
714 1752 : smearing = smear
715 : ELSE
716 : smearing = .FALSE.
717 : END IF
718 :
719 1752 : IF (PRESENT(kTS_sum)) THEN
720 1752 : entropic_term = kTS_sum
721 : ELSE
722 : entropic_term = 0.0_dp
723 : END IF
724 :
725 1752 : NULLIFY (energy)
726 1752 : CALL get_qs_env(qs_env, energy=energy)
727 1752 : CALL almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
728 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
729 1752 : print_active=.TRUE.)
730 :
731 : !! Add electronic entropy contribution if smearing is requested
732 : !! Previous QS entropy is replaced by the sum of the entropy for each spin
733 1752 : IF (smearing) THEN
734 20 : energy%total = energy%total - energy%kTS + entropic_term
735 : END IF
736 :
737 1752 : energy_total = energy%total
738 :
739 1752 : CALL timestop(handle)
740 :
741 1752 : END SUBROUTINE almo_dm_to_qs_ks
742 :
743 : ! **************************************************************************************************
744 : !> \brief uses the ALMO density matrix
745 : !> to compute ALMO KS matrix and the new energy
746 : !> \param qs_env ...
747 : !> \param matrix_p ...
748 : !> \param matrix_ks ...
749 : !> \param energy_total ...
750 : !> \param eps_filter ...
751 : !> \param mat_distr_aos ...
752 : !> \param smear ...
753 : !> \param kTS_sum ...
754 : !> \par History
755 : !> 2011.05 created [Rustam Z Khaliullin]
756 : !> 2018.09 smearing support [Ruben Staub]
757 : !> \author Rustam Z Khaliullin
758 : ! **************************************************************************************************
759 1752 : SUBROUTINE almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, &
760 : mat_distr_aos, smear, kTS_sum)
761 :
762 : TYPE(qs_environment_type), POINTER :: qs_env
763 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p, matrix_ks
764 : REAL(KIND=dp) :: energy_total, eps_filter
765 : INTEGER, INTENT(IN) :: mat_distr_aos
766 : LOGICAL, INTENT(IN), OPTIONAL :: smear
767 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS_sum
768 :
769 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_almo_ks'
770 :
771 : INTEGER :: handle, ispin, nspins
772 : LOGICAL :: smearing
773 : REAL(KIND=dp) :: entropic_term
774 1752 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_qs_ks
775 :
776 1752 : CALL timeset(routineN, handle)
777 :
778 1752 : IF (PRESENT(smear)) THEN
779 470 : smearing = smear
780 : ELSE
781 1282 : smearing = .FALSE.
782 : END IF
783 :
784 1752 : IF (PRESENT(kTS_sum)) THEN
785 470 : entropic_term = kTS_sum
786 : ELSE
787 1282 : entropic_term = 0.0_dp
788 : END IF
789 :
790 : ! update KS matrix in the QS env
791 : CALL almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, &
792 : smear=smearing, &
793 1752 : kTS_sum=entropic_term)
794 :
795 1752 : nspins = SIZE(matrix_ks)
796 :
797 : ! get KS matrix from the QS env and convert to the ALMO format
798 1752 : CALL get_qs_env(qs_env, matrix_ks=matrix_qs_ks)
799 3522 : DO ispin = 1, nspins
800 1770 : CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
801 3522 : CALL dbcsr_filter(matrix_ks(ispin), eps_filter)
802 : END DO
803 :
804 1752 : CALL timestop(handle)
805 :
806 1752 : END SUBROUTINE almo_dm_to_almo_ks
807 :
808 : ! **************************************************************************************************
809 : !> \brief update qs_env total energy
810 : !> \param qs_env ...
811 : !> \param energy ...
812 : !> \param energy_singles_corr ...
813 : !> \par History
814 : !> 2013.03 created [Rustam Z Khaliullin]
815 : !> \author Rustam Z Khaliullin
816 : ! **************************************************************************************************
817 112 : SUBROUTINE almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr)
818 :
819 : TYPE(qs_environment_type), POINTER :: qs_env
820 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: energy, energy_singles_corr
821 :
822 : TYPE(qs_energy_type), POINTER :: qs_energy
823 :
824 112 : CALL get_qs_env(qs_env, energy=qs_energy)
825 :
826 112 : IF (PRESENT(energy_singles_corr)) THEN
827 26 : qs_energy%singles_corr = energy_singles_corr
828 : ELSE
829 86 : qs_energy%singles_corr = 0.0_dp
830 : END IF
831 :
832 112 : IF (PRESENT(energy)) THEN
833 112 : qs_energy%total = energy
834 : END IF
835 :
836 112 : qs_energy%total = qs_energy%total + qs_energy%singles_corr
837 :
838 112 : END SUBROUTINE almo_scf_update_ks_energy
839 :
840 : ! **************************************************************************************************
841 : !> \brief Creates the matrix that imposes absolute locality on MOs
842 : !> \param qs_env ...
843 : !> \param almo_scf_env ...
844 : !> \par History
845 : !> 2011.11 created [Rustam Z. Khaliullin]
846 : !> \author Rustam Z. Khaliullin
847 : ! **************************************************************************************************
848 122 : SUBROUTINE almo_scf_construct_quencher(qs_env, almo_scf_env)
849 :
850 : TYPE(qs_environment_type), POINTER :: qs_env
851 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
852 :
853 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_construct_quencher'
854 :
855 : CHARACTER :: sym
856 : INTEGER :: col, contact_atom_1, contact_atom_2, domain_col, domain_map_local_entries, &
857 : domain_row, global_entries, global_list_length, grid1, GroupID, handle, hold, iatom, &
858 : iatom2, iblock_col, iblock_row, idomain, idomain2, ientry, igrid, ineig, ineighbor, &
859 : iNode, inode2, ipair, ispin, jatom, jatom2, jdomain2, local_list_length, &
860 : max_domain_neighbors, max_neig, mynode, nblkcols_tot, nblkrows_tot, nblks, ndomains, &
861 : neig_temp, nnode2, nNodes, row, unit_nr
862 122 : INTEGER, ALLOCATABLE, DIMENSION(:) :: current_number_neighbors, domain_entries_cpu, &
863 122 : domain_map_global, domain_map_local, first_atom_of_molecule, global_list, &
864 122 : last_atom_of_molecule, list_length_cpu, list_offset_cpu, local_list, offset_for_cpu
865 122 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: domain_grid, domain_neighbor_list, &
866 122 : domain_neighbor_list_excessive
867 122 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
868 : LOGICAL :: already_listed, block_active, &
869 : delayed_increment, found, &
870 : max_neig_fails, tr
871 : REAL(KIND=dp) :: contact1_radius, contact2_radius, &
872 : distance, distance_squared, overlap, &
873 : r0, r1, s0, s1, trial_distance_squared
874 122 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: new_block
875 : REAL(KIND=dp), DIMENSION(3) :: rab
876 122 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_old_block
877 : TYPE(cell_type), POINTER :: cell
878 : TYPE(cp_logger_type), POINTER :: logger
879 : TYPE(dbcsr_distribution_type) :: dist
880 122 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
881 : TYPE(dbcsr_type) :: matrix_s_sym
882 122 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
883 : TYPE(mp_comm_type) :: group
884 : TYPE(neighbor_list_iterator_p_type), &
885 122 : DIMENSION(:), POINTER :: nl_iterator, nl_iterator2
886 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
887 122 : POINTER :: sab_almo
888 122 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
889 :
890 122 : CALL timeset(routineN, handle)
891 :
892 : ! get a useful output_unit
893 122 : logger => cp_get_default_logger()
894 122 : IF (logger%para_env%is_source()) THEN
895 61 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
896 : ELSE
897 : unit_nr = -1
898 : END IF
899 :
900 122 : ndomains = almo_scf_env%ndomains
901 :
902 : CALL get_qs_env(qs_env=qs_env, &
903 : particle_set=particle_set, &
904 : molecule_set=molecule_set, &
905 : cell=cell, &
906 : matrix_s=matrix_s, &
907 122 : sab_almo=sab_almo)
908 :
909 : ! if we are dealing with molecules get info about them
910 122 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular .OR. &
911 : almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
912 366 : ALLOCATE (first_atom_of_molecule(almo_scf_env%nmolecules))
913 244 : ALLOCATE (last_atom_of_molecule(almo_scf_env%nmolecules))
914 : CALL get_molecule_set_info(molecule_set, &
915 : mol_to_first_atom=first_atom_of_molecule, &
916 122 : mol_to_last_atom=last_atom_of_molecule)
917 : END IF
918 :
919 : ! create a symmetrized copy of the ao overlap
920 : CALL dbcsr_create(matrix_s_sym, &
921 : template=almo_scf_env%matrix_s(1), &
922 122 : matrix_type=dbcsr_type_no_symmetry)
923 : CALL dbcsr_get_info(almo_scf_env%matrix_s(1), &
924 122 : matrix_type=sym)
925 122 : IF (sym == dbcsr_type_no_symmetry) THEN
926 0 : CALL dbcsr_copy(matrix_s_sym, almo_scf_env%matrix_s(1))
927 : ELSE
928 : CALL dbcsr_desymmetrize(almo_scf_env%matrix_s(1), &
929 122 : matrix_s_sym)
930 : END IF
931 :
932 494 : ALLOCATE (almo_scf_env%quench_t(almo_scf_env%nspins))
933 494 : ALLOCATE (almo_scf_env%domain_map(almo_scf_env%nspins))
934 :
935 250 : DO ispin = 1, almo_scf_env%nspins
936 :
937 : ! create the sparsity template for the occupied orbitals
938 : CALL matrix_almo_create(matrix_new=almo_scf_env%quench_t(ispin), &
939 : matrix_qs=matrix_s(1)%matrix, &
940 : almo_scf_env=almo_scf_env, &
941 : name_new="T_QUENCHER", &
942 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_occ], &
943 : symmetry_new=dbcsr_type_no_symmetry, &
944 : spin_key=ispin, &
945 128 : init_domains=.FALSE.)
946 :
947 : ! initialize distance quencher
948 128 : CALL dbcsr_work_create(almo_scf_env%quench_t(ispin), work_mutable=.TRUE.)
949 : CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist, &
950 128 : nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
951 128 : CALL dbcsr_distribution_get(dist, numnodes=nNodes, group=GroupID, mynode=mynode)
952 128 : CALL group%set_handle(groupid)
953 :
954 : ! create global atom neighbor list from the local lists
955 : ! first, calculate number of local pairs
956 128 : local_list_length = 0
957 128 : CALL neighbor_list_iterator_create(nl_iterator, sab_almo)
958 39611 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
959 : ! nnode - total number of neighbors for iatom
960 : ! inode - current neighbor count
961 : CALL get_iterator_info(nl_iterator, &
962 39483 : iatom=iatom2, jatom=jatom2, inode=inode2, nnode=nnode2)
963 39611 : IF (inode2 == 1) THEN
964 2049 : local_list_length = local_list_length + nnode2
965 : END IF
966 : END DO
967 128 : CALL neighbor_list_iterator_release(nl_iterator)
968 :
969 : ! second, extract the local list to an array
970 383 : ALLOCATE (local_list(2*local_list_length))
971 128 : local_list(:) = 0
972 128 : local_list_length = 0
973 128 : CALL neighbor_list_iterator_create(nl_iterator2, sab_almo)
974 39611 : DO WHILE (neighbor_list_iterate(nl_iterator2) == 0)
975 : CALL get_iterator_info(nl_iterator2, &
976 39483 : iatom=iatom2, jatom=jatom2)
977 39483 : local_list(2*local_list_length + 1) = iatom2
978 39483 : local_list(2*local_list_length + 2) = jatom2
979 39483 : local_list_length = local_list_length + 1
980 : END DO ! end loop over pairs of atoms
981 128 : CALL neighbor_list_iterator_release(nl_iterator2)
982 :
983 : ! third, communicate local length to the other nodes
984 512 : ALLOCATE (list_length_cpu(nNodes), list_offset_cpu(nNodes))
985 128 : CALL group%allgather(2*local_list_length, list_length_cpu)
986 :
987 : ! fourth, create a global list
988 128 : list_offset_cpu(1) = 0
989 256 : DO iNode = 2, nNodes
990 : list_offset_cpu(iNode) = list_offset_cpu(iNode - 1) + &
991 256 : list_length_cpu(iNode - 1)
992 : END DO
993 128 : global_list_length = list_offset_cpu(nNodes) + list_length_cpu(nNodes)
994 :
995 : ! fifth, communicate all list data
996 384 : ALLOCATE (global_list(global_list_length))
997 : CALL group%allgatherv(local_list, global_list, &
998 128 : list_length_cpu, list_offset_cpu)
999 128 : DEALLOCATE (list_length_cpu, list_offset_cpu)
1000 128 : DEALLOCATE (local_list)
1001 :
1002 : ! calculate maximum number of atoms surrounding the domain
1003 384 : ALLOCATE (current_number_neighbors(almo_scf_env%ndomains))
1004 128 : current_number_neighbors(:) = 0
1005 128 : global_list_length = global_list_length/2
1006 79094 : DO ipair = 1, global_list_length
1007 78966 : iatom2 = global_list(2*(ipair - 1) + 1)
1008 78966 : jatom2 = global_list(2*(ipair - 1) + 2)
1009 78966 : idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
1010 78966 : jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
1011 : ! add to the list
1012 78966 : current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
1013 : ! add j,i with i,j
1014 79094 : IF (idomain2 /= jdomain2) THEN
1015 63144 : current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
1016 : END IF
1017 : END DO
1018 962 : max_domain_neighbors = MAXVAL(current_number_neighbors)
1019 :
1020 : ! use the global atom neighbor list to create a global domain neighbor list
1021 512 : ALLOCATE (domain_neighbor_list_excessive(ndomains, max_domain_neighbors))
1022 962 : current_number_neighbors(:) = 1
1023 962 : DO ipair = 1, ndomains
1024 962 : domain_neighbor_list_excessive(ipair, 1) = ipair
1025 : END DO
1026 79094 : DO ipair = 1, global_list_length
1027 78966 : iatom2 = global_list(2*(ipair - 1) + 1)
1028 78966 : jatom2 = global_list(2*(ipair - 1) + 2)
1029 78966 : idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
1030 78966 : jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
1031 78966 : already_listed = .FALSE.
1032 325938 : DO ineighbor = 1, current_number_neighbors(idomain2)
1033 325938 : IF (domain_neighbor_list_excessive(idomain2, ineighbor) == jdomain2) THEN
1034 : already_listed = .TRUE.
1035 : EXIT
1036 : END IF
1037 : END DO
1038 79094 : IF (.NOT. already_listed) THEN
1039 : ! add to the list
1040 2722 : current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
1041 2722 : domain_neighbor_list_excessive(idomain2, current_number_neighbors(idomain2)) = jdomain2
1042 : ! add j,i with i,j
1043 2722 : IF (idomain2 /= jdomain2) THEN
1044 2722 : current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
1045 2722 : domain_neighbor_list_excessive(jdomain2, current_number_neighbors(jdomain2)) = idomain2
1046 : END IF
1047 : END IF
1048 : END DO ! end loop over pairs of atoms
1049 128 : DEALLOCATE (global_list)
1050 :
1051 962 : max_domain_neighbors = MAXVAL(current_number_neighbors)
1052 512 : ALLOCATE (domain_neighbor_list(ndomains, max_domain_neighbors))
1053 128 : domain_neighbor_list(:, :) = 0
1054 7248 : domain_neighbor_list(:, :) = domain_neighbor_list_excessive(:, 1:max_domain_neighbors)
1055 128 : DEALLOCATE (domain_neighbor_list_excessive)
1056 :
1057 384 : ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1058 384 : ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(max_domain_neighbors*ndomains, 2))
1059 12956 : almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1060 962 : almo_scf_env%domain_map(ispin)%index1(:) = 0
1061 128 : domain_map_local_entries = 0
1062 :
1063 : ! RZK-warning intermediate [0,1] quencher values are ill-defined
1064 : ! for molecules (not continuous and conceptually inadequate)
1065 :
1066 : CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), &
1067 128 : row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1068 : ! O(N) loop over domain pairs
1069 962 : DO row = 1, nblkrows_tot
1070 7240 : DO col = 1, current_number_neighbors(row)
1071 6278 : tr = .FALSE.
1072 6278 : iblock_row = row
1073 6278 : iblock_col = domain_neighbor_list(row, col)
1074 : CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin), &
1075 6278 : iblock_row, iblock_col, hold)
1076 :
1077 7112 : IF (hold == mynode) THEN
1078 :
1079 : ! Translate indices of distribution blocks to indices of domain blocks
1080 : ! Rows are AOs
1081 3139 : domain_row = almo_scf_env%domain_index_of_ao_block(iblock_row)
1082 : ! Columns are electrons (i.e. MOs)
1083 3139 : domain_col = almo_scf_env%domain_index_of_mo_block(iblock_col)
1084 :
1085 3139 : SELECT CASE (almo_scf_env%constraint_type)
1086 : CASE (almo_constraint_block_diagonal)
1087 :
1088 0 : block_active = .FALSE.
1089 : ! type of electron groups
1090 0 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1091 :
1092 : ! type of ao domains
1093 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1094 :
1095 : ! ao domains are molecular / electron groups are molecular
1096 0 : IF (domain_row == domain_col) THEN
1097 : block_active = .TRUE.
1098 : END IF
1099 :
1100 : ELSE ! ao domains are atomic
1101 :
1102 : ! ao domains are atomic / electron groups are molecular
1103 0 : CPABORT("Illegal: atomic domains and molecular groups")
1104 :
1105 : END IF
1106 :
1107 : ELSE ! electron groups are atomic
1108 :
1109 : ! type of ao domains
1110 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1111 :
1112 : ! ao domains are molecular / electron groups are atomic
1113 0 : CPABORT("Illegal: molecular domains and atomic groups")
1114 :
1115 : ELSE
1116 :
1117 : ! ao domains are atomic / electron groups are atomic
1118 0 : IF (domain_row == domain_col) THEN
1119 : block_active = .TRUE.
1120 : END IF
1121 :
1122 : END IF
1123 :
1124 : END IF ! end type of electron groups
1125 :
1126 0 : IF (block_active) THEN
1127 :
1128 0 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1129 0 : new_block(:, :) = 1.0_dp
1130 0 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1131 0 : DEALLOCATE (new_block)
1132 :
1133 0 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1134 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1135 : END IF
1136 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1137 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1138 0 : domain_map_local_entries = domain_map_local_entries + 1
1139 :
1140 : END IF
1141 :
1142 : CASE (almo_constraint_ao_overlap)
1143 :
1144 : ! type of electron groups
1145 0 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1146 :
1147 : ! type of ao domains
1148 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1149 :
1150 : ! ao domains are molecular / electron groups are molecular
1151 :
1152 : ! compute the maximum overlap between the atoms of the two molecules
1153 0 : CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1154 0 : IF (found) THEN
1155 0 : overlap = MAXVAL(ABS(p_old_block))
1156 : ELSE
1157 : overlap = 0.0_dp
1158 : END IF
1159 :
1160 : ELSE ! ao domains are atomic
1161 :
1162 : ! ao domains are atomic / electron groups are molecular
1163 : ! overlap_between_atom_and_molecule(atom=domain_row,molecule=domain_col)
1164 0 : CPABORT("atomic domains and molecular groups - NYI")
1165 :
1166 : END IF
1167 :
1168 : ELSE ! electron groups are atomic
1169 :
1170 : ! type of ao domains
1171 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1172 :
1173 : ! ao domains are molecular / electron groups are atomic
1174 : ! overlap_between_atom_and_molecule(atom=domain_col,molecule=domain_row)
1175 0 : CPABORT("molecular domains and atomic groups - NYI")
1176 :
1177 : ELSE
1178 :
1179 : ! ao domains are atomic / electron groups are atomic
1180 : ! compute max overlap between atoms: domain_row and domain_col
1181 0 : CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1182 0 : IF (found) THEN
1183 0 : overlap = MAXVAL(ABS(p_old_block))
1184 : ELSE
1185 : overlap = 0.0_dp
1186 : END IF
1187 :
1188 : END IF
1189 :
1190 : END IF ! end type of electron groups
1191 :
1192 0 : s0 = -LOG10(ABS(almo_scf_env%quencher_s0))
1193 0 : s1 = -LOG10(ABS(almo_scf_env%quencher_s1))
1194 0 : IF (overlap == 0.0_dp) THEN
1195 0 : overlap = -LOG10(ABS(almo_scf_env%eps_filter)) + 100.0_dp
1196 : ELSE
1197 0 : overlap = -LOG10(overlap)
1198 : END IF
1199 0 : IF (s0 < 0.0_dp) THEN
1200 0 : CPABORT("S0 is less than zero")
1201 : END IF
1202 0 : IF (s1 <= 0.0_dp) THEN
1203 0 : CPABORT("S1 is less than or equal to zero")
1204 : END IF
1205 0 : IF (s0 >= s1) THEN
1206 0 : CPABORT("S0 is greater than or equal to S1")
1207 : END IF
1208 :
1209 : ! Fill in non-zero blocks if AOs are close to the electron center
1210 0 : IF (overlap < s1) THEN
1211 0 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1212 0 : IF (overlap <= s0) THEN
1213 0 : new_block(:, :) = 1.0_dp
1214 : ELSE
1215 0 : new_block(:, :) = 1.0_dp/(1.0_dp + EXP(-(s0 - s1)/(s0 - overlap) - (s0 - s1)/(overlap - s1)))
1216 : END IF
1217 :
1218 0 : IF (ABS(new_block(1, 1)) > ABS(almo_scf_env%eps_filter)) THEN
1219 0 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1220 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1221 : END IF
1222 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1223 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1224 0 : domain_map_local_entries = domain_map_local_entries + 1
1225 : END IF
1226 :
1227 0 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1228 0 : DEALLOCATE (new_block)
1229 :
1230 : END IF
1231 :
1232 : CASE (almo_constraint_distance)
1233 :
1234 : ! type of electron groups
1235 3139 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1236 :
1237 : ! type of ao domains
1238 3139 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1239 :
1240 : ! ao domains are molecular / electron groups are molecular
1241 :
1242 : ! compute distance between molecules: domain_row and domain_col
1243 : ! distance between molecules is defined as the smallest
1244 : ! distance among all atom pairs
1245 3139 : IF (domain_row == domain_col) THEN
1246 417 : distance = 0.0_dp
1247 417 : contact_atom_1 = first_atom_of_molecule(domain_row)
1248 417 : contact_atom_2 = first_atom_of_molecule(domain_col)
1249 : ELSE
1250 2722 : distance_squared = 1.0E+100_dp
1251 2722 : contact_atom_1 = -1
1252 2722 : contact_atom_2 = -1
1253 9074 : DO iatom = first_atom_of_molecule(domain_row), last_atom_of_molecule(domain_row)
1254 26414 : DO jatom = first_atom_of_molecule(domain_col), last_atom_of_molecule(domain_col)
1255 17340 : rab(:) = pbc(particle_set(iatom)%r(:), particle_set(jatom)%r(:), cell)
1256 17340 : trial_distance_squared = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1257 23692 : IF (trial_distance_squared < distance_squared) THEN
1258 6391 : distance_squared = trial_distance_squared
1259 6391 : contact_atom_1 = iatom
1260 6391 : contact_atom_2 = jatom
1261 : END IF
1262 : END DO ! jatom
1263 : END DO ! iatom
1264 2722 : CPASSERT(contact_atom_1 > 0)
1265 2722 : distance = SQRT(distance_squared)
1266 : END IF
1267 :
1268 : ELSE ! ao domains are atomic
1269 :
1270 : ! ao domains are atomic / electron groups are molecular
1271 : !distance_between_atom_and_molecule(atom=domain_row,molecule=domain_col)
1272 0 : CPABORT("atomic domains and molecular groups - NYI")
1273 :
1274 : END IF
1275 :
1276 : ELSE ! electron groups are atomic
1277 :
1278 : ! type of ao domains
1279 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1280 :
1281 : ! ao domains are molecular / electron groups are atomic
1282 : !distance_between_atom_and_molecule(atom=domain_col,molecule=domain_row)
1283 0 : CPABORT("molecular domains and atomic groups - NYI")
1284 :
1285 : ELSE
1286 :
1287 : ! ao domains are atomic / electron groups are atomic
1288 : ! compute distance between atoms: domain_row and domain_col
1289 0 : rab(:) = pbc(particle_set(domain_row)%r(:), particle_set(domain_col)%r(:), cell)
1290 0 : distance = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1291 0 : contact_atom_1 = domain_row
1292 0 : contact_atom_2 = domain_col
1293 :
1294 : END IF
1295 :
1296 : END IF ! end type of electron groups
1297 :
1298 : ! get atomic radii to compute distance cutoff threshold
1299 3139 : IF (almo_scf_env%quencher_radius_type == do_bondparm_covalent) THEN
1300 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1301 0 : rcov=contact1_radius)
1302 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1303 0 : rcov=contact2_radius)
1304 3139 : ELSE IF (almo_scf_env%quencher_radius_type == do_bondparm_vdw) THEN
1305 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1306 3139 : rvdw=contact1_radius)
1307 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1308 3139 : rvdw=contact2_radius)
1309 : ELSE
1310 0 : CPABORT("Illegal quencher_radius_type")
1311 : END IF
1312 3139 : contact1_radius = cp_unit_to_cp2k(contact1_radius, "angstrom")
1313 3139 : contact2_radius = cp_unit_to_cp2k(contact2_radius, "angstrom")
1314 :
1315 : !RZK-warning the procedure is faulty for molecules:
1316 : ! the closest contacts should be found using
1317 : ! the element specific radii
1318 :
1319 : ! compute inner and outer cutoff radii
1320 3139 : r0 = almo_scf_env%quencher_r0_factor*(contact1_radius + contact2_radius)
1321 : !+almo_scf_env%quencher_r0_shift
1322 3139 : r1 = almo_scf_env%quencher_r1_factor*(contact1_radius + contact2_radius)
1323 : !+almo_scf_env%quencher_r1_shift
1324 :
1325 3139 : IF (r0 < 0.0_dp) THEN
1326 0 : CPABORT("R0 is less than zero")
1327 : END IF
1328 3139 : IF (r1 <= 0.0_dp) THEN
1329 0 : CPABORT("R1 is less than or equal to zero")
1330 : END IF
1331 3139 : IF (r0 > r1) THEN
1332 0 : CPABORT("R0 is greater than or equal to R1")
1333 : END IF
1334 :
1335 : ! Fill in non-zero blocks if AOs are close to the electron center
1336 3139 : IF (distance < r1) THEN
1337 8740 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1338 2185 : IF (distance <= r0) THEN
1339 101919 : new_block(:, :) = 1.0_dp
1340 : ELSE
1341 : ! remove the intermediate values from the quencher temporarily
1342 0 : CPABORT("")
1343 0 : new_block(:, :) = 1.0_dp/(1.0_dp + EXP((r1 - r0)/(r0 - distance) + (r1 - r0)/(r1 - distance)))
1344 : END IF
1345 :
1346 2185 : IF (ABS(new_block(1, 1)) > ABS(almo_scf_env%eps_filter)) THEN
1347 2185 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1348 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1349 : END IF
1350 2185 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1351 2185 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1352 2185 : domain_map_local_entries = domain_map_local_entries + 1
1353 : END IF
1354 :
1355 2185 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1356 2185 : DEALLOCATE (new_block)
1357 : END IF
1358 :
1359 : CASE DEFAULT
1360 3139 : CPABORT("Illegal constraint type")
1361 : END SELECT
1362 :
1363 : END IF ! mynode
1364 :
1365 : END DO
1366 : END DO ! end O(N) loop over pairs
1367 :
1368 128 : DEALLOCATE (domain_neighbor_list)
1369 128 : DEALLOCATE (current_number_neighbors)
1370 :
1371 128 : CALL dbcsr_finalize(almo_scf_env%quench_t(ispin))
1372 :
1373 : CALL dbcsr_filter(almo_scf_env%quench_t(ispin), &
1374 128 : almo_scf_env%eps_filter)
1375 :
1376 : ! check that both domain_map and quench_t have the same number of entries
1377 128 : nblks = dbcsr_get_num_blocks(almo_scf_env%quench_t(ispin))
1378 128 : IF (nblks /= domain_map_local_entries) THEN
1379 0 : CPABORT("number of blocks is wrong")
1380 : END IF
1381 :
1382 : ! first, communicate map sizes on the other nodes
1383 384 : ALLOCATE (domain_entries_cpu(nNodes), offset_for_cpu(nNodes))
1384 128 : CALL group%allgather(2*domain_map_local_entries, domain_entries_cpu)
1385 :
1386 : ! second, create
1387 128 : offset_for_cpu(1) = 0
1388 256 : DO iNode = 2, nNodes
1389 : offset_for_cpu(iNode) = offset_for_cpu(iNode - 1) + &
1390 256 : domain_entries_cpu(iNode - 1)
1391 : END DO
1392 128 : global_entries = offset_for_cpu(nNodes) + domain_entries_cpu(nNodes)
1393 :
1394 : ! communicate all entries
1395 384 : ALLOCATE (domain_map_global(global_entries))
1396 508 : ALLOCATE (domain_map_local(2*domain_map_local_entries))
1397 2313 : DO ientry = 1, domain_map_local_entries
1398 2185 : domain_map_local(2*(ientry - 1) + 1) = almo_scf_env%domain_map(ispin)%pairs(ientry, 1)
1399 2313 : domain_map_local(2*ientry) = almo_scf_env%domain_map(ispin)%pairs(ientry, 2)
1400 : END DO
1401 : CALL group%allgatherv(domain_map_local, domain_map_global, &
1402 128 : domain_entries_cpu, offset_for_cpu)
1403 128 : DEALLOCATE (domain_entries_cpu, offset_for_cpu)
1404 128 : DEALLOCATE (domain_map_local)
1405 :
1406 128 : DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
1407 128 : DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
1408 256 : ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1409 512 : ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(global_entries/2, 2))
1410 9124 : almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1411 962 : almo_scf_env%domain_map(ispin)%index1(:) = 0
1412 :
1413 : ! unpack the received data into a local variable
1414 : ! since we do not know the maximum global number of neighbors
1415 : ! try one. if fails increase the maximum number and try again
1416 : ! until it succeeds
1417 : max_neig = max_domain_neighbors
1418 : max_neig_fails = .TRUE.
1419 256 : max_neig_loop: DO WHILE (max_neig_fails)
1420 512 : ALLOCATE (domain_grid(almo_scf_env%ndomains, 0:max_neig))
1421 128 : domain_grid(:, :) = 0
1422 : ! init the number of collected neighbors
1423 962 : domain_grid(:, 0) = 1
1424 : ! loop over the records
1425 128 : global_entries = global_entries/2
1426 4498 : DO ientry = 1, global_entries
1427 : ! get the center
1428 4370 : grid1 = domain_map_global(2*ientry)
1429 : ! get the neighbor
1430 4370 : ineig = domain_map_global(2*(ientry - 1) + 1)
1431 : ! check boundaries
1432 4370 : IF (domain_grid(grid1, 0) > max_neig) THEN
1433 : ! this neighbor will overstep the boundaries
1434 : ! stop the trial and increase the max number of neighbors
1435 0 : DEALLOCATE (domain_grid)
1436 0 : max_neig = max_neig*2
1437 0 : CYCLE max_neig_loop
1438 : END IF
1439 : ! for the current center loop over the collected neighbors
1440 : ! to insert the current record in a numerical order
1441 : delayed_increment = .FALSE.
1442 19016 : DO igrid = 1, domain_grid(grid1, 0)
1443 : ! compare the current neighbor with that already in the 'book'
1444 19016 : IF (ineig < domain_grid(grid1, igrid)) THEN
1445 : ! if this one is smaller then insert it here and pick up the one
1446 : ! from the book to continue inserting
1447 3180 : neig_temp = ineig
1448 3180 : ineig = domain_grid(grid1, igrid)
1449 3180 : domain_grid(grid1, igrid) = neig_temp
1450 : ELSE
1451 11466 : IF (domain_grid(grid1, igrid) == 0) THEN
1452 : ! got the empty slot now - insert the record
1453 4370 : domain_grid(grid1, igrid) = ineig
1454 : ! increase the record counter but do it outside the loop
1455 4370 : delayed_increment = .TRUE.
1456 : END IF
1457 : END IF
1458 : END DO
1459 4498 : IF (delayed_increment) THEN
1460 4370 : domain_grid(grid1, 0) = domain_grid(grid1, 0) + 1
1461 : ELSE
1462 : ! should not be here - all records must be inserted
1463 0 : CPABORT("all records must be inserted")
1464 : END IF
1465 : END DO
1466 128 : max_neig_fails = .FALSE.
1467 : END DO max_neig_loop
1468 128 : DEALLOCATE (domain_map_global)
1469 :
1470 128 : ientry = 1
1471 962 : DO idomain = 1, almo_scf_env%ndomains
1472 5204 : DO ineig = 1, domain_grid(idomain, 0) - 1
1473 4370 : almo_scf_env%domain_map(ispin)%pairs(ientry, 1) = domain_grid(idomain, ineig)
1474 4370 : almo_scf_env%domain_map(ispin)%pairs(ientry, 2) = idomain
1475 5204 : ientry = ientry + 1
1476 : END DO
1477 962 : almo_scf_env%domain_map(ispin)%index1(idomain) = ientry
1478 : END DO
1479 378 : DEALLOCATE (domain_grid)
1480 :
1481 : END DO ! ispin
1482 122 : IF (almo_scf_env%nspins == 2) THEN
1483 : CALL dbcsr_copy(almo_scf_env%quench_t(2), &
1484 6 : almo_scf_env%quench_t(1))
1485 : almo_scf_env%domain_map(2)%pairs(:, :) = &
1486 50 : almo_scf_env%domain_map(1)%pairs(:, :)
1487 : almo_scf_env%domain_map(2)%index1(:) = &
1488 18 : almo_scf_env%domain_map(1)%index1(:)
1489 : END IF
1490 :
1491 122 : CALL dbcsr_release(matrix_s_sym)
1492 :
1493 122 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular .OR. &
1494 : almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1495 122 : DEALLOCATE (first_atom_of_molecule)
1496 122 : DEALLOCATE (last_atom_of_molecule)
1497 : END IF
1498 :
1499 : !mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(&
1500 : ! dbcsr_distribution(almo_scf_env%quench_t(ispin))))
1501 : !CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist, &
1502 : ! nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
1503 : !DO row = 1, nblkrows_tot
1504 : ! DO col = 1, nblkcols_tot
1505 : ! tr = .FALSE.
1506 : ! iblock_row = row
1507 : ! iblock_col = col
1508 : ! CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin),&
1509 : ! iblock_row, iblock_col, tr, hold)
1510 : ! CALL dbcsr_get_block_p(almo_scf_env%quench_t(ispin),&
1511 : ! row, col, p_old_block, found)
1512 : ! write(*,*) "RST_NOTE:", mynode, row, col, hold, found
1513 : ! ENDDO
1514 : !ENDDO
1515 :
1516 122 : CALL timestop(handle)
1517 :
1518 244 : END SUBROUTINE almo_scf_construct_quencher
1519 :
1520 : ! *****************************************************************************
1521 : !> \brief Compute matrix W (energy-weighted density matrix) that is needed
1522 : !> for the evaluation of forces
1523 : !> \param matrix_w ...
1524 : !> \param almo_scf_env ...
1525 : !> \par History
1526 : !> 2015.03 created [Rustam Z. Khaliullin]
1527 : !> \author Rustam Z. Khaliullin
1528 : ! **************************************************************************************************
1529 66 : SUBROUTINE calculate_w_matrix_almo(matrix_w, almo_scf_env)
1530 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
1531 : TYPE(almo_scf_env_type) :: almo_scf_env
1532 :
1533 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_almo'
1534 :
1535 : INTEGER :: handle, ispin
1536 : REAL(KIND=dp) :: scaling
1537 : TYPE(dbcsr_type) :: tmp_nn1, tmp_no1, tmp_oo1, tmp_oo2
1538 :
1539 66 : CALL timeset(routineN, handle)
1540 :
1541 66 : IF (almo_scf_env%nspins == 1) THEN
1542 66 : scaling = 2.0_dp
1543 : ELSE
1544 0 : scaling = 1.0_dp
1545 : END IF
1546 :
1547 132 : DO ispin = 1, almo_scf_env%nspins
1548 :
1549 : CALL dbcsr_create(tmp_nn1, template=almo_scf_env%matrix_s(1), &
1550 66 : matrix_type=dbcsr_type_no_symmetry)
1551 : CALL dbcsr_create(tmp_no1, template=almo_scf_env%matrix_t(ispin), &
1552 66 : matrix_type=dbcsr_type_no_symmetry)
1553 : CALL dbcsr_create(tmp_oo1, template=almo_scf_env%matrix_sigma_inv(ispin), &
1554 66 : matrix_type=dbcsr_type_no_symmetry)
1555 : CALL dbcsr_create(tmp_oo2, template=almo_scf_env%matrix_sigma_inv(ispin), &
1556 66 : matrix_type=dbcsr_type_no_symmetry)
1557 :
1558 66 : CALL dbcsr_copy(tmp_nn1, almo_scf_env%matrix_ks(ispin))
1559 : ! 1. TMP_NO1=F.T
1560 : CALL dbcsr_multiply("N", "N", scaling, tmp_nn1, almo_scf_env%matrix_t(ispin), &
1561 66 : 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1562 : ! 2. TMP_OO1=T^(tr).TMP_NO1=T^(tr).(FT)
1563 : CALL dbcsr_multiply("T", "N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_no1, &
1564 66 : 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1565 : ! 3. TMP_OO2=TMP_OO1.siginv=(T^(tr)FT).siginv
1566 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_oo1, almo_scf_env%matrix_sigma_inv(ispin), &
1567 66 : 0.0_dp, tmp_oo2, filter_eps=almo_scf_env%eps_filter)
1568 : ! 4. TMP_OO1=siginv.TMP_OO2=siginv.(T^(tr)FTsiginv)
1569 : CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_sigma_inv(ispin), tmp_oo2, &
1570 66 : 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1571 : ! 5. TMP_NO1=T.TMP_OO1.=T.(siginvT^(tr)FTsiginv)
1572 : CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_oo1, &
1573 66 : 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1574 : ! 6. TMP_NN1=TMP_NO1.T^(tr)=(TsiginvT^(tr)FTsiginv).T^(tr)=RFR
1575 : CALL dbcsr_multiply("N", "T", 1.0_dp, tmp_no1, almo_scf_env%matrix_t(ispin), &
1576 66 : 0.0_dp, tmp_nn1, filter_eps=almo_scf_env%eps_filter)
1577 66 : CALL matrix_almo_to_qs(tmp_nn1, matrix_w(ispin)%matrix, almo_scf_env%mat_distr_aos)
1578 :
1579 66 : CALL dbcsr_release(tmp_nn1)
1580 66 : CALL dbcsr_release(tmp_no1)
1581 66 : CALL dbcsr_release(tmp_oo1)
1582 132 : CALL dbcsr_release(tmp_oo2)
1583 :
1584 : END DO
1585 :
1586 66 : CALL timestop(handle)
1587 :
1588 66 : END SUBROUTINE calculate_w_matrix_almo
1589 :
1590 : END MODULE almo_scf_qs
1591 :
|
