Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 3284 : 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 3284 : INTEGER, DIMENSION(:), POINTER :: blk_distr, blk_sizes, block_sizes_new, col_blk_size, &
132 3284 : 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 3284 : 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 3284 : 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 3284 : CALL dbcsr_get_info(matrix_qs, distribution=dist_qs)
163 :
164 3284 : natoms = almo_scf_env%natoms
165 3284 : nmols = almo_scf_env%nmolecules
166 :
167 9852 : DO dimen = 1, 2 ! 1 - row, 2 - column dimension
168 :
169 : ! distribution pattern is the same for all matrix types (ao, occ, virt)
170 6568 : IF (dimen == 1) THEN !rows
171 3284 : CALL dbcsr_distribution_get(dist_qs, row_dist=blk_distr)
172 : ELSE !columns
173 3284 : CALL dbcsr_distribution_get(dist_qs, col_dist=blk_distr)
174 : END IF
175 :
176 6568 : 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 2508 : CALL dbcsr_get_info(matrix_qs, row_blk_size=blk_sizes)
180 :
181 : ! atomic clustering of AOs
182 2508 : 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 2508 : ELSE IF (almo_scf_env%mat_distr_aos == almo_mat_distr_molecular) THEN
188 12540 : ALLOCATE (block_sizes_new(nmols), distr_new_array(nmols))
189 20130 : block_sizes_new(:) = 0
190 45936 : 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 45936 : blk_sizes(iatom)
194 : END DO
195 20130 : DO imol = 1, nmols
196 : distr_new_array(imol) = &
197 20130 : 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 4060 : 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 4060 : 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 4060 : ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN
224 4060 : nlength = nmols
225 12180 : ALLOCATE (block_sizes_new(nlength))
226 4060 : IF (size_keys(dimen) == almo_mat_dim_occ) THEN
227 18520 : 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 18520 : WHERE (block_sizes_new == 0) block_sizes_new = 1
231 1740 : 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 1740 : ELSE IF (size_keys(dimen) == almo_mat_dim_virt_full) THEN
234 5556 : block_sizes_new(:) = almo_scf_env%nvirt_full_of_domain(:, spin_key)
235 1044 : ELSE IF (size_keys(dimen) == almo_mat_dim_virt) THEN
236 8334 : 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 12180 : ALLOCATE (distr_new_array(nlength))
250 : ! atomic clustering
251 4060 : IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN
252 0 : distr_new_array(:) = blk_distr(:)
253 : ! molecular clustering
254 4060 : ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN
255 32410 : DO imol = 1, nmols
256 : distr_new_array(imol) = &
257 32410 : 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 9852 : IF (dimen == 1) THEN !rows
264 3284 : row_sizes_new => block_sizes_new
265 3284 : row_distr_new => distr_new_array
266 : ELSE !columns
267 3284 : col_sizes_new => block_sizes_new
268 3284 : 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 3284 : reuse_arrays=.TRUE.)
276 :
277 : ! Create the matrix
278 : CALL dbcsr_create(matrix_new, name_new, &
279 : dist_new, symmetry_new, &
280 3284 : row_sizes_new, col_sizes_new, reuse_arrays=.TRUE.)
281 3284 : 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 3284 : IF (init_domains) THEN
286 :
287 1312 : CALL dbcsr_distribution_get(dist_new, mynode=mynode)
288 1312 : CALL dbcsr_work_create(matrix_new, work_mutable=.TRUE.)
289 : CALL dbcsr_get_info(matrix_new, nblkrows_total=nblkrows_tot, &
290 1312 : 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 10528 : DO row = 1, nblkrows_tot
351 9216 : tr = .FALSE.
352 9216 : iblock_row = row
353 9216 : iblock_col = row
354 9216 : CALL dbcsr_get_stored_coordinates(matrix_new, iblock_row, iblock_col, hold)
355 :
356 10528 : IF (hold == mynode) THEN
357 :
358 13824 : active = .TRUE.
359 :
360 : one_dim_is_mo = .FALSE.
361 13824 : DO dimen = 1, 2 ! 1 - row, 2 - column dimension
362 13824 : IF (size_keys(dimen) == almo_mat_dim_occ) one_dim_is_mo = .TRUE.
363 : END DO
364 4608 : IF (one_dim_is_mo) THEN
365 1620 : IF (almo_scf_env%nocc_of_domain(row, spin_key) == 0) active = .FALSE.
366 : END IF
367 :
368 4608 : one_dim_is_mo = .FALSE.
369 13824 : DO dimen = 1, 2
370 13824 : IF (size_keys(dimen) == almo_mat_dim_virt) one_dim_is_mo = .TRUE.
371 : END DO
372 4608 : IF (one_dim_is_mo) THEN
373 405 : IF (almo_scf_env%nvirt_of_domain(row, spin_key) == 0) active = .FALSE.
374 : END IF
375 :
376 4608 : one_dim_is_mo = .FALSE.
377 13824 : DO dimen = 1, 2
378 13824 : IF (size_keys(dimen) == almo_mat_dim_virt_disc) one_dim_is_mo = .TRUE.
379 : END DO
380 4608 : 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 4608 : one_dim_is_mo = .FALSE.
385 13824 : DO dimen = 1, 2
386 13824 : IF (size_keys(dimen) == almo_mat_dim_virt_full) one_dim_is_mo = .TRUE.
387 : END DO
388 4608 : IF (one_dim_is_mo) THEN
389 405 : IF (almo_scf_env%nvirt_full_of_domain(row, spin_key) == 0) active = .FALSE.
390 : END IF
391 :
392 4574 : IF (active) THEN
393 17136 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
394 837014 : new_block(:, :) = 1.0_dp
395 4284 : CALL dbcsr_put_block(matrix_new, iblock_row, iblock_col, new_block)
396 4284 : 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 3284 : CALL dbcsr_finalize(matrix_new)
406 :
407 3284 : CALL timestop(handle)
408 :
409 6568 : 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 2026 : 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 2026 : CALL timeset(routineN, handle)
431 : !RZK-warning if it's not a N(AO)xN(AO) matrix then stop
432 :
433 2026 : 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 2026 : CALL dbcsr_create(matrix_qs_nosym, template=matrix_qs, matrix_type=dbcsr_type_no_symmetry)
440 2026 : 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 2026 : CALL dbcsr_set(matrix_almo, 0.0_dp)
450 2026 : CALL dbcsr_complete_redistribute(matrix_qs_nosym, matrix_almo)
451 2026 : CALL dbcsr_release(matrix_qs_nosym)
452 :
453 : CASE DEFAULT
454 2026 : CPABORT("")
455 : END SELECT
456 :
457 2026 : CALL timestop(handle)
458 :
459 2026 : 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 1826 : 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 1826 : CALL timeset(routineN, handle)
480 : ! RZK-warning if it's not a N(AO)xN(AO) matrix then stop
481 :
482 1826 : 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 1826 : CALL dbcsr_create(matrix_almo_redist, template=matrix_qs)
487 1826 : CALL dbcsr_complete_redistribute(matrix_almo, matrix_almo_redist)
488 1826 : CALL dbcsr_set(matrix_qs, 0.0_dp)
489 1826 : CALL dbcsr_copy(matrix_qs, matrix_almo_redist, keep_sparsity=.TRUE.)
490 1826 : CALL dbcsr_release(matrix_almo_redist)
491 : CASE DEFAULT
492 1826 : CPABORT("")
493 : END SELECT
494 :
495 1826 : CALL timestop(handle)
496 :
497 1826 : 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 116 : 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 116 : 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 116 : POINTER :: sab_orb
523 : TYPE(qs_ks_env_type), POINTER :: ks_env
524 :
525 116 : CALL timeset(routineN, handle)
526 :
527 116 : 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 116 : sab_orb=sab_orb)
536 :
537 116 : nspin = dft_control%nspins
538 :
539 : ! create matrix_ks in the QS env if necessary
540 116 : 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 232 : DO ispin = 1, nspin
554 116 : CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
555 232 : CALL dbcsr_filter(matrix_ks(ispin), eps_filter)
556 : END DO
557 :
558 116 : CALL timestop(handle)
559 :
560 116 : 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 348 : 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 116 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
582 : TYPE(qs_scf_env_type), POINTER :: scf_env
583 :
584 116 : CALL timeset(routineN, handle)
585 :
586 : ! create and init scf_env (this is necessary to return MOs to qs)
587 116 : NULLIFY (mos, fm_struct_tmp, scf_env)
588 116 : ALLOCATE (scf_env)
589 116 : CALL scf_env_create(scf_env)
590 :
591 : !CALL qs_scf_env_initialize(qs_env, scf_env)
592 116 : CALL set_qs_env(qs_env, scf_env=scf_env)
593 116 : CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
594 :
595 116 : 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 232 : DO ispin = 1, almo_scf_env%nspins
599 :
600 : ! Currently only fm version of mo_set is usable.
601 : ! First transform the matrix_t to fm version
602 : ! Empty the containers to prevent memory leaks
603 116 : CALL deallocate_mo_set(mos(ispin))
604 : CALL allocate_mo_set(mo_set=mos(ispin), &
605 : nao=nrow_fm, &
606 : nmo=ncol_fm, &
607 : nelectron=almo_scf_env%nelectrons_total, &
608 : n_el_f=REAL(almo_scf_env%nelectrons_total, dp), &
609 : maxocc=2.0_dp, &
610 116 : flexible_electron_count=dft_control%relax_multiplicity)
611 :
612 : CALL cp_fm_struct_create(fm_struct_tmp, nrow_global=nrow_fm, ncol_global=ncol_fm, &
613 : context=almo_scf_env%blacs_env, &
614 116 : para_env=almo_scf_env%para_env)
615 :
616 116 : CALL cp_fm_create(mo_fm_copy, fm_struct_tmp, name="t_orthogonal_converted_to_fm")
617 116 : CALL cp_fm_struct_release(fm_struct_tmp)
618 : !CALL copy_dbcsr_to_fm(almo_scf_env%matrix_t(ispin), mo_fm_copy)
619 :
620 116 : CALL init_mo_set(mos(ispin), fm_ref=mo_fm_copy, name='fm_mo')
621 :
622 348 : CALL cp_fm_release(mo_fm_copy)
623 :
624 : END DO
625 :
626 116 : CALL timestop(handle)
627 :
628 116 : END SUBROUTINE construct_qs_mos
629 :
630 : ! **************************************************************************************************
631 : !> \brief return density matrix to the qs_env
632 : !> \param qs_env ...
633 : !> \param matrix_p ...
634 : !> \param mat_distr_aos ...
635 : !> \par History
636 : !> 2011.05 created [Rustam Z Khaliullin]
637 : !> \author Rustam Z Khaliullin
638 : ! **************************************************************************************************
639 1760 : SUBROUTINE almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
640 : TYPE(qs_environment_type), POINTER :: qs_env
641 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
642 : INTEGER, INTENT(IN) :: mat_distr_aos
643 :
644 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_qs_env'
645 :
646 : INTEGER :: handle, ispin, nspins
647 1760 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
648 : TYPE(qs_rho_type), POINTER :: rho
649 :
650 1760 : CALL timeset(routineN, handle)
651 :
652 1760 : NULLIFY (rho, rho_ao)
653 1760 : nspins = SIZE(matrix_p)
654 1760 : CALL get_qs_env(qs_env, rho=rho)
655 1760 : CALL qs_rho_get(rho, rho_ao=rho_ao)
656 :
657 : ! set the new density matrix
658 3520 : DO ispin = 1, nspins
659 : CALL matrix_almo_to_qs(matrix_p(ispin), &
660 : rho_ao(ispin)%matrix, &
661 3520 : mat_distr_aos)
662 : END DO
663 1760 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
664 1760 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE.)
665 :
666 1760 : CALL timestop(handle)
667 :
668 1760 : END SUBROUTINE almo_dm_to_qs_env
669 :
670 : ! **************************************************************************************************
671 : !> \brief uses the ALMO density matrix
672 : !> to compute KS matrix (inside QS environment) and the new energy
673 : !> \param qs_env ...
674 : !> \param matrix_p ...
675 : !> \param energy_total ...
676 : !> \param mat_distr_aos ...
677 : !> \param smear ...
678 : !> \param kTS_sum ...
679 : !> \par History
680 : !> 2011.05 created [Rustam Z Khaliullin]
681 : !> 2018.09 smearing support [Ruben Staub]
682 : !> \author Rustam Z Khaliullin
683 : ! **************************************************************************************************
684 1734 : SUBROUTINE almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, smear, kTS_sum)
685 : TYPE(qs_environment_type), POINTER :: qs_env
686 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p
687 : REAL(KIND=dp) :: energy_total
688 : INTEGER, INTENT(IN) :: mat_distr_aos
689 : LOGICAL, INTENT(IN), OPTIONAL :: smear
690 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS_sum
691 :
692 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_qs_ks'
693 :
694 : INTEGER :: handle
695 : LOGICAL :: smearing
696 : REAL(KIND=dp) :: entropic_term
697 : TYPE(qs_energy_type), POINTER :: energy
698 :
699 1734 : CALL timeset(routineN, handle)
700 :
701 1734 : IF (PRESENT(smear)) THEN
702 1734 : smearing = smear
703 : ELSE
704 : smearing = .FALSE.
705 : END IF
706 :
707 1734 : IF (PRESENT(kTS_sum)) THEN
708 1734 : entropic_term = kTS_sum
709 : ELSE
710 : entropic_term = 0.0_dp
711 : END IF
712 :
713 1734 : NULLIFY (energy)
714 1734 : CALL get_qs_env(qs_env, energy=energy)
715 1734 : CALL almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
716 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.FALSE., &
717 1734 : print_active=.TRUE.)
718 :
719 : !! Add electronic entropy contribution if smearing is requested
720 : !! Previous QS entropy is replaced by the sum of the entropy for each spin
721 1734 : IF (smearing) THEN
722 20 : energy%total = energy%total - energy%kTS + entropic_term
723 : END IF
724 :
725 1734 : energy_total = energy%total
726 :
727 1734 : CALL timestop(handle)
728 :
729 1734 : END SUBROUTINE almo_dm_to_qs_ks
730 :
731 : ! **************************************************************************************************
732 : !> \brief uses the ALMO density matrix
733 : !> to compute ALMO KS matrix and the new energy
734 : !> \param qs_env ...
735 : !> \param matrix_p ...
736 : !> \param matrix_ks ...
737 : !> \param energy_total ...
738 : !> \param eps_filter ...
739 : !> \param mat_distr_aos ...
740 : !> \param smear ...
741 : !> \param kTS_sum ...
742 : !> \par History
743 : !> 2011.05 created [Rustam Z Khaliullin]
744 : !> 2018.09 smearing support [Ruben Staub]
745 : !> \author Rustam Z Khaliullin
746 : ! **************************************************************************************************
747 1734 : SUBROUTINE almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, &
748 : mat_distr_aos, smear, kTS_sum)
749 :
750 : TYPE(qs_environment_type), POINTER :: qs_env
751 : TYPE(dbcsr_type), DIMENSION(:) :: matrix_p, matrix_ks
752 : REAL(KIND=dp) :: energy_total, eps_filter
753 : INTEGER, INTENT(IN) :: mat_distr_aos
754 : LOGICAL, INTENT(IN), OPTIONAL :: smear
755 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kTS_sum
756 :
757 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_dm_to_almo_ks'
758 :
759 : INTEGER :: handle, ispin, nspins
760 : LOGICAL :: smearing
761 : REAL(KIND=dp) :: entropic_term
762 1734 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_qs_ks
763 :
764 1734 : CALL timeset(routineN, handle)
765 :
766 1734 : IF (PRESENT(smear)) THEN
767 464 : smearing = smear
768 : ELSE
769 1270 : smearing = .FALSE.
770 : END IF
771 :
772 1734 : IF (PRESENT(kTS_sum)) THEN
773 464 : entropic_term = kTS_sum
774 : ELSE
775 1270 : entropic_term = 0.0_dp
776 : END IF
777 :
778 : ! update KS matrix in the QS env
779 : CALL almo_dm_to_qs_ks(qs_env, matrix_p, energy_total, mat_distr_aos, &
780 : smear=smearing, &
781 1734 : kTS_sum=entropic_term)
782 :
783 1734 : nspins = SIZE(matrix_ks)
784 :
785 : ! get KS matrix from the QS env and convert to the ALMO format
786 1734 : CALL get_qs_env(qs_env, matrix_ks=matrix_qs_ks)
787 3468 : DO ispin = 1, nspins
788 1734 : CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, matrix_ks(ispin), mat_distr_aos)
789 3468 : CALL dbcsr_filter(matrix_ks(ispin), eps_filter)
790 : END DO
791 :
792 1734 : CALL timestop(handle)
793 :
794 1734 : END SUBROUTINE almo_dm_to_almo_ks
795 :
796 : ! **************************************************************************************************
797 : !> \brief update qs_env total energy
798 : !> \param qs_env ...
799 : !> \param energy ...
800 : !> \param energy_singles_corr ...
801 : !> \par History
802 : !> 2013.03 created [Rustam Z Khaliullin]
803 : !> \author Rustam Z Khaliullin
804 : ! **************************************************************************************************
805 106 : SUBROUTINE almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr)
806 :
807 : TYPE(qs_environment_type), POINTER :: qs_env
808 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: energy, energy_singles_corr
809 :
810 : TYPE(qs_energy_type), POINTER :: qs_energy
811 :
812 106 : CALL get_qs_env(qs_env, energy=qs_energy)
813 :
814 106 : IF (PRESENT(energy_singles_corr)) THEN
815 26 : qs_energy%singles_corr = energy_singles_corr
816 : ELSE
817 80 : qs_energy%singles_corr = 0.0_dp
818 : END IF
819 :
820 106 : IF (PRESENT(energy)) THEN
821 106 : qs_energy%total = energy
822 : END IF
823 :
824 106 : qs_energy%total = qs_energy%total + qs_energy%singles_corr
825 :
826 106 : END SUBROUTINE almo_scf_update_ks_energy
827 :
828 : ! **************************************************************************************************
829 : !> \brief Creates the matrix that imposes absolute locality on MOs
830 : !> \param qs_env ...
831 : !> \param almo_scf_env ...
832 : !> \par History
833 : !> 2011.11 created [Rustam Z. Khaliullin]
834 : !> \author Rustam Z. Khaliullin
835 : ! **************************************************************************************************
836 116 : SUBROUTINE almo_scf_construct_quencher(qs_env, almo_scf_env)
837 :
838 : TYPE(qs_environment_type), POINTER :: qs_env
839 : TYPE(almo_scf_env_type), INTENT(INOUT) :: almo_scf_env
840 :
841 : CHARACTER(len=*), PARAMETER :: routineN = 'almo_scf_construct_quencher'
842 :
843 : CHARACTER :: sym
844 : INTEGER :: col, contact_atom_1, contact_atom_2, domain_col, domain_map_local_entries, &
845 : domain_row, global_entries, global_list_length, grid1, GroupID, handle, hold, iatom, &
846 : iatom2, iblock_col, iblock_row, idomain, idomain2, ientry, igrid, ineig, ineighbor, &
847 : iNode, inode2, ipair, ispin, jatom, jatom2, jdomain2, local_list_length, &
848 : max_domain_neighbors, max_neig, mynode, nblkcols_tot, nblkrows_tot, nblks, ndomains, &
849 : neig_temp, nnode2, nNodes, row, unit_nr
850 116 : INTEGER, ALLOCATABLE, DIMENSION(:) :: current_number_neighbors, domain_entries_cpu, &
851 116 : domain_map_global, domain_map_local, first_atom_of_molecule, global_list, &
852 116 : last_atom_of_molecule, list_length_cpu, list_offset_cpu, local_list, offset_for_cpu
853 116 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: domain_grid, domain_neighbor_list, &
854 116 : domain_neighbor_list_excessive
855 116 : INTEGER, DIMENSION(:), POINTER :: col_blk_size, row_blk_size
856 : LOGICAL :: already_listed, block_active, &
857 : delayed_increment, found, &
858 : max_neig_fails, tr
859 : REAL(KIND=dp) :: contact1_radius, contact2_radius, &
860 : distance, distance_squared, overlap, &
861 : r0, r1, s0, s1, trial_distance_squared
862 116 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: new_block
863 : REAL(KIND=dp), DIMENSION(3) :: rab
864 116 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_old_block
865 : TYPE(cell_type), POINTER :: cell
866 : TYPE(cp_logger_type), POINTER :: logger
867 : TYPE(dbcsr_distribution_type) :: dist
868 116 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
869 : TYPE(dbcsr_type) :: matrix_s_sym
870 116 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
871 : TYPE(mp_comm_type) :: group
872 : TYPE(neighbor_list_iterator_p_type), &
873 116 : DIMENSION(:), POINTER :: nl_iterator, nl_iterator2
874 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
875 116 : POINTER :: sab_almo
876 116 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
877 :
878 116 : CALL timeset(routineN, handle)
879 :
880 : ! get a useful output_unit
881 116 : logger => cp_get_default_logger()
882 116 : IF (logger%para_env%is_source()) THEN
883 58 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
884 : ELSE
885 : unit_nr = -1
886 : END IF
887 :
888 116 : ndomains = almo_scf_env%ndomains
889 :
890 : CALL get_qs_env(qs_env=qs_env, &
891 : particle_set=particle_set, &
892 : molecule_set=molecule_set, &
893 : cell=cell, &
894 : matrix_s=matrix_s, &
895 116 : sab_almo=sab_almo)
896 :
897 : ! if we are dealing with molecules get info about them
898 116 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular .OR. &
899 : almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
900 348 : ALLOCATE (first_atom_of_molecule(almo_scf_env%nmolecules))
901 232 : ALLOCATE (last_atom_of_molecule(almo_scf_env%nmolecules))
902 : CALL get_molecule_set_info(molecule_set, &
903 : mol_to_first_atom=first_atom_of_molecule, &
904 116 : mol_to_last_atom=last_atom_of_molecule)
905 : END IF
906 :
907 : ! create a symmetrized copy of the ao overlap
908 : CALL dbcsr_create(matrix_s_sym, &
909 : template=almo_scf_env%matrix_s(1), &
910 116 : matrix_type=dbcsr_type_no_symmetry)
911 : CALL dbcsr_get_info(almo_scf_env%matrix_s(1), &
912 116 : matrix_type=sym)
913 116 : IF (sym == dbcsr_type_no_symmetry) THEN
914 0 : CALL dbcsr_copy(matrix_s_sym, almo_scf_env%matrix_s(1))
915 : ELSE
916 : CALL dbcsr_desymmetrize(almo_scf_env%matrix_s(1), &
917 116 : matrix_s_sym)
918 : END IF
919 :
920 464 : ALLOCATE (almo_scf_env%quench_t(almo_scf_env%nspins))
921 464 : ALLOCATE (almo_scf_env%domain_map(almo_scf_env%nspins))
922 :
923 : !DO ispin=1,almo_scf_env%nspins
924 116 : ispin = 1
925 :
926 : ! create the sparsity template for the occupied orbitals
927 : CALL matrix_almo_create(matrix_new=almo_scf_env%quench_t(ispin), &
928 : matrix_qs=matrix_s(1)%matrix, &
929 : almo_scf_env=almo_scf_env, &
930 : name_new="T_QUENCHER", &
931 : size_keys=[almo_mat_dim_aobasis, almo_mat_dim_occ], &
932 : symmetry_new=dbcsr_type_no_symmetry, &
933 : spin_key=ispin, &
934 116 : init_domains=.FALSE.)
935 :
936 : ! initialize distance quencher
937 116 : CALL dbcsr_work_create(almo_scf_env%quench_t(ispin), work_mutable=.TRUE.)
938 : CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist, &
939 116 : nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
940 116 : CALL dbcsr_distribution_get(dist, numnodes=nNodes, group=GroupID, mynode=mynode)
941 116 : CALL group%set_handle(groupid)
942 :
943 : ! create global atom neighbor list from the local lists
944 : ! first, calculate number of local pairs
945 116 : local_list_length = 0
946 116 : CALL neighbor_list_iterator_create(nl_iterator, sab_almo)
947 39355 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
948 : ! nnode - total number of neighbors for iatom
949 : ! inode - current neighbor count
950 : CALL get_iterator_info(nl_iterator, &
951 39239 : iatom=iatom2, jatom=jatom2, inode=inode2, nnode=nnode2)
952 : !WRITE(*,*) "GET INFO: ",iatom2, jatom2, inode2, nnode2
953 39355 : IF (inode2 == 1) THEN
954 2001 : local_list_length = local_list_length + nnode2
955 : END IF
956 : END DO
957 116 : CALL neighbor_list_iterator_release(nl_iterator)
958 :
959 : ! second, extract the local list to an array
960 347 : ALLOCATE (local_list(2*local_list_length))
961 78594 : local_list(:) = 0
962 116 : local_list_length = 0
963 116 : CALL neighbor_list_iterator_create(nl_iterator2, sab_almo)
964 39355 : DO WHILE (neighbor_list_iterate(nl_iterator2) == 0)
965 : CALL get_iterator_info(nl_iterator2, &
966 39239 : iatom=iatom2, jatom=jatom2)
967 39239 : local_list(2*local_list_length + 1) = iatom2
968 39239 : local_list(2*local_list_length + 2) = jatom2
969 39239 : local_list_length = local_list_length + 1
970 : END DO ! end loop over pairs of atoms
971 116 : CALL neighbor_list_iterator_release(nl_iterator2)
972 :
973 : ! third, communicate local length to the other nodes
974 464 : ALLOCATE (list_length_cpu(nNodes), list_offset_cpu(nNodes))
975 116 : CALL group%allgather(2*local_list_length, list_length_cpu)
976 :
977 : ! fourth, create a global list
978 116 : list_offset_cpu(1) = 0
979 232 : DO iNode = 2, nNodes
980 : list_offset_cpu(iNode) = list_offset_cpu(iNode - 1) + &
981 232 : list_length_cpu(iNode - 1)
982 : END DO
983 116 : global_list_length = list_offset_cpu(nNodes) + list_length_cpu(nNodes)
984 :
985 : ! fifth, communicate all list data
986 348 : ALLOCATE (global_list(global_list_length))
987 : CALL group%allgatherv(local_list, global_list, &
988 116 : list_length_cpu, list_offset_cpu)
989 116 : DEALLOCATE (list_length_cpu, list_offset_cpu)
990 116 : DEALLOCATE (local_list)
991 :
992 : ! calculate maximum number of atoms surrounding the domain
993 348 : ALLOCATE (current_number_neighbors(almo_scf_env%ndomains))
994 926 : current_number_neighbors(:) = 0
995 116 : global_list_length = global_list_length/2
996 78594 : DO ipair = 1, global_list_length
997 78478 : iatom2 = global_list(2*(ipair - 1) + 1)
998 78478 : jatom2 = global_list(2*(ipair - 1) + 2)
999 78478 : idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
1000 78478 : jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
1001 : ! add to the list
1002 78478 : current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
1003 : ! add j,i with i,j
1004 78594 : IF (idomain2 /= jdomain2) THEN
1005 62940 : current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
1006 : END IF
1007 : END DO
1008 926 : max_domain_neighbors = MAXVAL(current_number_neighbors)
1009 :
1010 : ! use the global atom neighbor list to create a global domain neighbor list
1011 464 : ALLOCATE (domain_neighbor_list_excessive(ndomains, max_domain_neighbors))
1012 926 : current_number_neighbors(:) = 1
1013 926 : DO ipair = 1, ndomains
1014 926 : domain_neighbor_list_excessive(ipair, 1) = ipair
1015 : END DO
1016 78594 : DO ipair = 1, global_list_length
1017 78478 : iatom2 = global_list(2*(ipair - 1) + 1)
1018 78478 : jatom2 = global_list(2*(ipair - 1) + 2)
1019 78478 : idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
1020 78478 : jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
1021 78478 : already_listed = .FALSE.
1022 325246 : DO ineighbor = 1, current_number_neighbors(idomain2)
1023 325246 : IF (domain_neighbor_list_excessive(idomain2, ineighbor) == jdomain2) THEN
1024 : already_listed = .TRUE.
1025 : EXIT
1026 : END IF
1027 : END DO
1028 78594 : IF (.NOT. already_listed) THEN
1029 : ! add to the list
1030 2710 : current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
1031 2710 : domain_neighbor_list_excessive(idomain2, current_number_neighbors(idomain2)) = jdomain2
1032 : ! add j,i with i,j
1033 2710 : IF (idomain2 /= jdomain2) THEN
1034 2710 : current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
1035 2710 : domain_neighbor_list_excessive(jdomain2, current_number_neighbors(jdomain2)) = idomain2
1036 : END IF
1037 : END IF
1038 : END DO ! end loop over pairs of atoms
1039 116 : DEALLOCATE (global_list)
1040 :
1041 926 : max_domain_neighbors = MAXVAL(current_number_neighbors)
1042 464 : ALLOCATE (domain_neighbor_list(ndomains, max_domain_neighbors))
1043 7164 : domain_neighbor_list(:, :) = 0
1044 7164 : domain_neighbor_list(:, :) = domain_neighbor_list_excessive(:, 1:max_domain_neighbors)
1045 116 : DEALLOCATE (domain_neighbor_list_excessive)
1046 :
1047 348 : ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1048 348 : ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(max_domain_neighbors*ndomains, 2))
1049 12824 : almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1050 926 : almo_scf_env%domain_map(ispin)%index1(:) = 0
1051 116 : domain_map_local_entries = 0
1052 :
1053 : ! RZK-warning intermediate [0,1] quencher values are ill-defined
1054 : ! for molecules (not continuous and conceptually inadequate)
1055 :
1056 : CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), &
1057 116 : row_blk_size=row_blk_size, col_blk_size=col_blk_size)
1058 : ! O(N) loop over domain pairs
1059 926 : DO row = 1, nblkrows_tot
1060 7156 : DO col = 1, current_number_neighbors(row)
1061 6230 : tr = .FALSE.
1062 6230 : iblock_row = row
1063 6230 : iblock_col = domain_neighbor_list(row, col)
1064 : CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin), &
1065 6230 : iblock_row, iblock_col, hold)
1066 :
1067 7040 : IF (hold == mynode) THEN
1068 :
1069 : ! Translate indices of distribution blocks to indices of domain blocks
1070 : ! Rows are AOs
1071 3115 : domain_row = almo_scf_env%domain_index_of_ao_block(iblock_row)
1072 : ! Columns are electrons (i.e. MOs)
1073 3115 : domain_col = almo_scf_env%domain_index_of_mo_block(iblock_col)
1074 :
1075 3115 : SELECT CASE (almo_scf_env%constraint_type)
1076 : CASE (almo_constraint_block_diagonal)
1077 :
1078 0 : block_active = .FALSE.
1079 : ! type of electron groups
1080 0 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1081 :
1082 : ! type of ao domains
1083 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1084 :
1085 : ! ao domains are molecular / electron groups are molecular
1086 0 : IF (domain_row == domain_col) THEN
1087 : block_active = .TRUE.
1088 : END IF
1089 :
1090 : ELSE ! ao domains are atomic
1091 :
1092 : ! ao domains are atomic / electron groups are molecular
1093 0 : CPABORT("Illegal: atomic domains and molecular groups")
1094 :
1095 : END IF
1096 :
1097 : ELSE ! electron groups are atomic
1098 :
1099 : ! type of ao domains
1100 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1101 :
1102 : ! ao domains are molecular / electron groups are atomic
1103 0 : CPABORT("Illegal: molecular domains and atomic groups")
1104 :
1105 : ELSE
1106 :
1107 : ! ao domains are atomic / electron groups are atomic
1108 0 : IF (domain_row == domain_col) THEN
1109 : block_active = .TRUE.
1110 : END IF
1111 :
1112 : END IF
1113 :
1114 : END IF ! end type of electron groups
1115 :
1116 0 : IF (block_active) THEN
1117 :
1118 0 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1119 0 : new_block(:, :) = 1.0_dp
1120 0 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1121 0 : DEALLOCATE (new_block)
1122 :
1123 0 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1124 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1125 : END IF
1126 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1127 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1128 0 : domain_map_local_entries = domain_map_local_entries + 1
1129 :
1130 : END IF
1131 :
1132 : CASE (almo_constraint_ao_overlap)
1133 :
1134 : ! type of electron groups
1135 0 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1136 :
1137 : ! type of ao domains
1138 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1139 :
1140 : ! ao domains are molecular / electron groups are molecular
1141 :
1142 : ! compute the maximum overlap between the atoms of the two molecules
1143 0 : CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1144 0 : IF (found) THEN
1145 : ! CPErrorMessage(cp_failure_level,routineP,"S block not found")
1146 : ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
1147 0 : overlap = MAXVAL(ABS(p_old_block))
1148 : ELSE
1149 : overlap = 0.0_dp
1150 : END IF
1151 :
1152 : ELSE ! ao domains are atomic
1153 :
1154 : ! ao domains are atomic / electron groups are molecular
1155 : ! overlap_between_atom_and_molecule(atom=domain_row,molecule=domain_col)
1156 0 : CPABORT("atomic domains and molecular groups - NYI")
1157 :
1158 : END IF
1159 :
1160 : ELSE ! electron groups are atomic
1161 :
1162 : ! type of ao domains
1163 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1164 :
1165 : ! ao domains are molecular / electron groups are atomic
1166 : ! overlap_between_atom_and_molecule(atom=domain_col,molecule=domain_row)
1167 0 : CPABORT("molecular domains and atomic groups - NYI")
1168 :
1169 : ELSE
1170 :
1171 : ! ao domains are atomic / electron groups are atomic
1172 : ! compute max overlap between atoms: domain_row and domain_col
1173 0 : CALL dbcsr_get_block_p(matrix_s_sym, iblock_row, iblock_col, p_old_block, found)
1174 0 : IF (found) THEN
1175 : ! CPErrorMessage(cp_failure_level,routineP,"S block not found")
1176 : ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
1177 0 : overlap = MAXVAL(ABS(p_old_block))
1178 : ELSE
1179 : overlap = 0.0_dp
1180 : END IF
1181 :
1182 : END IF
1183 :
1184 : END IF ! end type of electron groups
1185 :
1186 0 : s0 = -LOG10(ABS(almo_scf_env%quencher_s0))
1187 0 : s1 = -LOG10(ABS(almo_scf_env%quencher_s1))
1188 0 : IF (overlap == 0.0_dp) THEN
1189 0 : overlap = -LOG10(ABS(almo_scf_env%eps_filter)) + 100.0_dp
1190 : ELSE
1191 0 : overlap = -LOG10(overlap)
1192 : END IF
1193 0 : IF (s0 < 0.0_dp) THEN
1194 0 : CPABORT("S0 is less than zero")
1195 : END IF
1196 0 : IF (s1 <= 0.0_dp) THEN
1197 0 : CPABORT("S1 is less than or equal to zero")
1198 : END IF
1199 0 : IF (s0 >= s1) THEN
1200 0 : CPABORT("S0 is greater than or equal to S1")
1201 : END IF
1202 :
1203 : ! Fill in non-zero blocks if AOs are close to the electron center
1204 0 : IF (overlap < s1) THEN
1205 0 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1206 0 : IF (overlap <= s0) THEN
1207 0 : new_block(:, :) = 1.0_dp
1208 : !WRITE(*,'(A15,2I7,3F8.3,E11.3)') "INTRA-BLOCKS: ",&
1209 : ! iblock_col, iblock_row, s0, s1, overlap, new_block(1,1)
1210 : ELSE
1211 0 : new_block(:, :) = 1.0_dp/(1.0_dp + EXP(-(s0 - s1)/(s0 - overlap) - (s0 - s1)/(overlap - s1)))
1212 : !WRITE(*,'(A15,2I7,3F8.3,E11.3)') "INTER-BLOCKS: ",&
1213 : ! iblock_col, iblock_row, s0, s1, overlap, new_block(1,1)
1214 : END IF
1215 :
1216 0 : IF (ABS(new_block(1, 1)) > ABS(almo_scf_env%eps_filter)) THEN
1217 0 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1218 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1219 : END IF
1220 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1221 0 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1222 0 : domain_map_local_entries = domain_map_local_entries + 1
1223 : END IF
1224 :
1225 0 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1226 0 : DEALLOCATE (new_block)
1227 :
1228 : END IF
1229 :
1230 : CASE (almo_constraint_distance)
1231 :
1232 : ! type of electron groups
1233 3115 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1234 :
1235 : ! type of ao domains
1236 3115 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1237 :
1238 : ! ao domains are molecular / electron groups are molecular
1239 :
1240 : ! compute distance between molecules: domain_row and domain_col
1241 : ! distance between molecules is defined as the smallest
1242 : ! distance among all atom pairs
1243 3115 : IF (domain_row == domain_col) THEN
1244 405 : distance = 0.0_dp
1245 405 : contact_atom_1 = first_atom_of_molecule(domain_row)
1246 405 : contact_atom_2 = first_atom_of_molecule(domain_col)
1247 : ELSE
1248 2710 : distance_squared = 1.0E+100_dp
1249 2710 : contact_atom_1 = -1
1250 2710 : contact_atom_2 = -1
1251 9034 : DO iatom = first_atom_of_molecule(domain_row), last_atom_of_molecule(domain_row)
1252 26310 : DO jatom = first_atom_of_molecule(domain_col), last_atom_of_molecule(domain_col)
1253 17276 : rab(:) = pbc(particle_set(iatom)%r(:), particle_set(jatom)%r(:), cell)
1254 17276 : trial_distance_squared = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1255 23600 : IF (trial_distance_squared < distance_squared) THEN
1256 6363 : distance_squared = trial_distance_squared
1257 6363 : contact_atom_1 = iatom
1258 6363 : contact_atom_2 = jatom
1259 : END IF
1260 : END DO ! jatom
1261 : END DO ! iatom
1262 2710 : CPASSERT(contact_atom_1 > 0)
1263 2710 : distance = SQRT(distance_squared)
1264 : END IF
1265 :
1266 : ELSE ! ao domains are atomic
1267 :
1268 : ! ao domains are atomic / electron groups are molecular
1269 : !distance_between_atom_and_molecule(atom=domain_row,molecule=domain_col)
1270 0 : CPABORT("atomic domains and molecular groups - NYI")
1271 :
1272 : END IF
1273 :
1274 : ELSE ! electron groups are atomic
1275 :
1276 : ! type of ao domains
1277 0 : IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1278 :
1279 : ! ao domains are molecular / electron groups are atomic
1280 : !distance_between_atom_and_molecule(atom=domain_col,molecule=domain_row)
1281 0 : CPABORT("molecular domains and atomic groups - NYI")
1282 :
1283 : ELSE
1284 :
1285 : ! ao domains are atomic / electron groups are atomic
1286 : ! compute distance between atoms: domain_row and domain_col
1287 0 : rab(:) = pbc(particle_set(domain_row)%r(:), particle_set(domain_col)%r(:), cell)
1288 0 : distance = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1289 0 : contact_atom_1 = domain_row
1290 0 : contact_atom_2 = domain_col
1291 :
1292 : END IF
1293 :
1294 : END IF ! end type of electron groups
1295 :
1296 : ! get atomic radii to compute distance cutoff threshold
1297 3115 : IF (almo_scf_env%quencher_radius_type == do_bondparm_covalent) THEN
1298 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1299 0 : rcov=contact1_radius)
1300 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1301 0 : rcov=contact2_radius)
1302 3115 : ELSE IF (almo_scf_env%quencher_radius_type == do_bondparm_vdw) THEN
1303 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_1)%atomic_kind, &
1304 3115 : rvdw=contact1_radius)
1305 : CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1306 3115 : rvdw=contact2_radius)
1307 : ELSE
1308 0 : CPABORT("Illegal quencher_radius_type")
1309 : END IF
1310 3115 : contact1_radius = cp_unit_to_cp2k(contact1_radius, "angstrom")
1311 3115 : contact2_radius = cp_unit_to_cp2k(contact2_radius, "angstrom")
1312 :
1313 : !RZK-warning the procedure is faulty for molecules:
1314 : ! the closest contacts should be found using
1315 : ! the element specific radii
1316 :
1317 : ! compute inner and outer cutoff radii
1318 3115 : r0 = almo_scf_env%quencher_r0_factor*(contact1_radius + contact2_radius)
1319 : !+almo_scf_env%quencher_r0_shift
1320 3115 : r1 = almo_scf_env%quencher_r1_factor*(contact1_radius + contact2_radius)
1321 : !+almo_scf_env%quencher_r1_shift
1322 :
1323 3115 : IF (r0 < 0.0_dp) THEN
1324 0 : CPABORT("R0 is less than zero")
1325 : END IF
1326 3115 : IF (r1 <= 0.0_dp) THEN
1327 0 : CPABORT("R1 is less than or equal to zero")
1328 : END IF
1329 3115 : IF (r0 > r1) THEN
1330 0 : CPABORT("R0 is greater than or equal to R1")
1331 : END IF
1332 :
1333 : ! Fill in non-zero blocks if AOs are close to the electron center
1334 3115 : IF (distance < r1) THEN
1335 8676 : ALLOCATE (new_block(row_blk_size(iblock_row), col_blk_size(iblock_col)))
1336 2169 : IF (distance <= r0) THEN
1337 101401 : new_block(:, :) = 1.0_dp
1338 : !WRITE(*,'(A15,2I7,5F8.3,E11.3)') "INTRA-BLOCKS: ",&
1339 : ! iblock_col, iblock_row, contact1_radius,&
1340 : ! contact2_radius, r0, r1, distance, new_block(1,1)
1341 : ELSE
1342 : ! remove the intermediate values from the quencher temporarily
1343 0 : CPABORT("")
1344 0 : new_block(:, :) = 1.0_dp/(1.0_dp + EXP((r1 - r0)/(r0 - distance) + (r1 - r0)/(r1 - distance)))
1345 : !WRITE(*,'(A15,2I7,5F8.3,E11.3)') "INTER-BLOCKS: ",&
1346 : ! iblock_col, iblock_row, contact1_radius,&
1347 : ! contact2_radius, r0, r1, distance, new_block(1,1)
1348 : END IF
1349 :
1350 2169 : IF (ABS(new_block(1, 1)) > ABS(almo_scf_env%eps_filter)) THEN
1351 2169 : IF (domain_map_local_entries >= max_domain_neighbors*almo_scf_env%ndomains) THEN
1352 0 : CPABORT("weird... max_domain_neighbors is exceeded")
1353 : END IF
1354 2169 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1355 2169 : almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1356 2169 : domain_map_local_entries = domain_map_local_entries + 1
1357 : END IF
1358 :
1359 2169 : CALL dbcsr_put_block(almo_scf_env%quench_t(ispin), iblock_row, iblock_col, new_block)
1360 2169 : DEALLOCATE (new_block)
1361 : END IF
1362 :
1363 : CASE DEFAULT
1364 3115 : CPABORT("Illegal constraint type")
1365 : END SELECT
1366 :
1367 : END IF ! mynode
1368 :
1369 : END DO
1370 : END DO ! end O(N) loop over pairs
1371 :
1372 116 : DEALLOCATE (domain_neighbor_list)
1373 116 : DEALLOCATE (current_number_neighbors)
1374 :
1375 116 : CALL dbcsr_finalize(almo_scf_env%quench_t(ispin))
1376 : !CALL dbcsr_scale(almo_scf_env%quench_t(ispin),&
1377 : ! almo_scf_env%envelope_amplitude)
1378 : CALL dbcsr_filter(almo_scf_env%quench_t(ispin), &
1379 116 : almo_scf_env%eps_filter)
1380 :
1381 : ! check that both domain_map and quench_t have the same number of entries
1382 116 : nblks = dbcsr_get_num_blocks(almo_scf_env%quench_t(ispin))
1383 116 : IF (nblks /= domain_map_local_entries) THEN
1384 0 : CPABORT("number of blocks is wrong")
1385 : END IF
1386 :
1387 : ! first, communicate map sizes on the other nodes
1388 348 : ALLOCATE (domain_entries_cpu(nNodes), offset_for_cpu(nNodes))
1389 116 : CALL group%allgather(2*domain_map_local_entries, domain_entries_cpu)
1390 :
1391 : ! second, create
1392 116 : offset_for_cpu(1) = 0
1393 232 : DO iNode = 2, nNodes
1394 : offset_for_cpu(iNode) = offset_for_cpu(iNode - 1) + &
1395 232 : domain_entries_cpu(iNode - 1)
1396 : END DO
1397 116 : global_entries = offset_for_cpu(nNodes) + domain_entries_cpu(nNodes)
1398 :
1399 : ! communicate all entries
1400 348 : ALLOCATE (domain_map_global(global_entries))
1401 460 : ALLOCATE (domain_map_local(2*domain_map_local_entries))
1402 2285 : DO ientry = 1, domain_map_local_entries
1403 2169 : domain_map_local(2*(ientry - 1) + 1) = almo_scf_env%domain_map(ispin)%pairs(ientry, 1)
1404 2285 : domain_map_local(2*ientry) = almo_scf_env%domain_map(ispin)%pairs(ientry, 2)
1405 : END DO
1406 : CALL group%allgatherv(domain_map_local, domain_map_global, &
1407 116 : domain_entries_cpu, offset_for_cpu)
1408 116 : DEALLOCATE (domain_entries_cpu, offset_for_cpu)
1409 116 : DEALLOCATE (domain_map_local)
1410 :
1411 116 : DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
1412 116 : DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
1413 232 : ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1414 464 : ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(global_entries/2, 2))
1415 9024 : almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1416 926 : almo_scf_env%domain_map(ispin)%index1(:) = 0
1417 :
1418 : ! unpack the received data into a local variable
1419 : ! since we do not know the maximum global number of neighbors
1420 : ! try one. if fails increase the maximum number and try again
1421 : ! until it succeeds
1422 : max_neig = max_domain_neighbors
1423 : max_neig_fails = .TRUE.
1424 232 : max_neig_loop: DO WHILE (max_neig_fails)
1425 464 : ALLOCATE (domain_grid(almo_scf_env%ndomains, 0:max_neig))
1426 8090 : domain_grid(:, :) = 0
1427 : ! init the number of collected neighbors
1428 926 : domain_grid(:, 0) = 1
1429 : ! loop over the records
1430 116 : global_entries = global_entries/2
1431 4454 : DO ientry = 1, global_entries
1432 : ! get the center
1433 4338 : grid1 = domain_map_global(2*ientry)
1434 : ! get the neighbor
1435 4338 : ineig = domain_map_global(2*(ientry - 1) + 1)
1436 : ! check boundaries
1437 4338 : IF (domain_grid(grid1, 0) > max_neig) THEN
1438 : ! this neighbor will overstep the boundaries
1439 : ! stop the trial and increase the max number of neighbors
1440 0 : DEALLOCATE (domain_grid)
1441 0 : max_neig = max_neig*2
1442 0 : CYCLE max_neig_loop
1443 : END IF
1444 : ! for the current center loop over the collected neighbors
1445 : ! to insert the current record in a numerical order
1446 : delayed_increment = .FALSE.
1447 18944 : DO igrid = 1, domain_grid(grid1, 0)
1448 : ! compare the current neighbor with that already in the 'book'
1449 18944 : IF (ineig < domain_grid(grid1, igrid)) THEN
1450 : ! if this one is smaller then insert it here and pick up the one
1451 : ! from the book to continue inserting
1452 3172 : neig_temp = ineig
1453 3172 : ineig = domain_grid(grid1, igrid)
1454 3172 : domain_grid(grid1, igrid) = neig_temp
1455 : ELSE
1456 11434 : IF (domain_grid(grid1, igrid) == 0) THEN
1457 : ! got the empty slot now - insert the record
1458 4338 : domain_grid(grid1, igrid) = ineig
1459 : ! increase the record counter but do it outside the loop
1460 4338 : delayed_increment = .TRUE.
1461 : END IF
1462 : END IF
1463 : END DO
1464 4454 : IF (delayed_increment) THEN
1465 4338 : domain_grid(grid1, 0) = domain_grid(grid1, 0) + 1
1466 : ELSE
1467 : ! should not be here - all records must be inserted
1468 0 : CPABORT("all records must be inserted")
1469 : END IF
1470 : END DO
1471 116 : max_neig_fails = .FALSE.
1472 : END DO max_neig_loop
1473 116 : DEALLOCATE (domain_map_global)
1474 :
1475 116 : ientry = 1
1476 926 : DO idomain = 1, almo_scf_env%ndomains
1477 5148 : DO ineig = 1, domain_grid(idomain, 0) - 1
1478 4338 : almo_scf_env%domain_map(ispin)%pairs(ientry, 1) = domain_grid(idomain, ineig)
1479 4338 : almo_scf_env%domain_map(ispin)%pairs(ientry, 2) = idomain
1480 5148 : ientry = ientry + 1
1481 : END DO
1482 926 : almo_scf_env%domain_map(ispin)%index1(idomain) = ientry
1483 : END DO
1484 116 : DEALLOCATE (domain_grid)
1485 :
1486 : !ENDDO ! ispin
1487 116 : IF (almo_scf_env%nspins == 2) THEN
1488 : CALL dbcsr_copy(almo_scf_env%quench_t(2), &
1489 0 : almo_scf_env%quench_t(1))
1490 : almo_scf_env%domain_map(2)%pairs(:, :) = &
1491 0 : almo_scf_env%domain_map(1)%pairs(:, :)
1492 : almo_scf_env%domain_map(2)%index1(:) = &
1493 0 : almo_scf_env%domain_map(1)%index1(:)
1494 : END IF
1495 :
1496 116 : CALL dbcsr_release(matrix_s_sym)
1497 :
1498 116 : IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular .OR. &
1499 : almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1500 116 : DEALLOCATE (first_atom_of_molecule)
1501 116 : DEALLOCATE (last_atom_of_molecule)
1502 : END IF
1503 :
1504 : !mynode = dbcsr_mp_mynode(dbcsr_distribution_mp(&
1505 : ! dbcsr_distribution(almo_scf_env%quench_t(ispin))))
1506 : !CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist, &
1507 : ! nblkrows_total=nblkrows_tot, nblkcols_total=nblkcols_tot)
1508 : !DO row = 1, nblkrows_tot
1509 : ! DO col = 1, nblkcols_tot
1510 : ! tr = .FALSE.
1511 : ! iblock_row = row
1512 : ! iblock_col = col
1513 : ! CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin),&
1514 : ! iblock_row, iblock_col, tr, hold)
1515 : ! CALL dbcsr_get_block_p(almo_scf_env%quench_t(ispin),&
1516 : ! row, col, p_old_block, found)
1517 : ! write(*,*) "RST_NOTE:", mynode, row, col, hold, found
1518 : ! ENDDO
1519 : !ENDDO
1520 :
1521 116 : CALL timestop(handle)
1522 :
1523 348 : END SUBROUTINE almo_scf_construct_quencher
1524 :
1525 : ! *****************************************************************************
1526 : !> \brief Compute matrix W (energy-weighted density matrix) that is needed
1527 : !> for the evaluation of forces
1528 : !> \param matrix_w ...
1529 : !> \param almo_scf_env ...
1530 : !> \par History
1531 : !> 2015.03 created [Rustam Z. Khaliullin]
1532 : !> \author Rustam Z. Khaliullin
1533 : ! **************************************************************************************************
1534 66 : SUBROUTINE calculate_w_matrix_almo(matrix_w, almo_scf_env)
1535 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_w
1536 : TYPE(almo_scf_env_type) :: almo_scf_env
1537 :
1538 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_w_matrix_almo'
1539 :
1540 : INTEGER :: handle, ispin
1541 : REAL(KIND=dp) :: scaling
1542 : TYPE(dbcsr_type) :: tmp_nn1, tmp_no1, tmp_oo1, tmp_oo2
1543 :
1544 66 : CALL timeset(routineN, handle)
1545 :
1546 66 : IF (almo_scf_env%nspins == 1) THEN
1547 66 : scaling = 2.0_dp
1548 : ELSE
1549 0 : scaling = 1.0_dp
1550 : END IF
1551 :
1552 132 : DO ispin = 1, almo_scf_env%nspins
1553 :
1554 : CALL dbcsr_create(tmp_nn1, template=almo_scf_env%matrix_s(1), &
1555 66 : matrix_type=dbcsr_type_no_symmetry)
1556 : CALL dbcsr_create(tmp_no1, template=almo_scf_env%matrix_t(ispin), &
1557 66 : matrix_type=dbcsr_type_no_symmetry)
1558 : CALL dbcsr_create(tmp_oo1, template=almo_scf_env%matrix_sigma_inv(ispin), &
1559 66 : matrix_type=dbcsr_type_no_symmetry)
1560 : CALL dbcsr_create(tmp_oo2, template=almo_scf_env%matrix_sigma_inv(ispin), &
1561 66 : matrix_type=dbcsr_type_no_symmetry)
1562 :
1563 66 : CALL dbcsr_copy(tmp_nn1, almo_scf_env%matrix_ks(ispin))
1564 : ! 1. TMP_NO1=F.T
1565 : CALL dbcsr_multiply("N", "N", scaling, tmp_nn1, almo_scf_env%matrix_t(ispin), &
1566 66 : 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1567 : ! 2. TMP_OO1=T^(tr).TMP_NO1=T^(tr).(FT)
1568 : CALL dbcsr_multiply("T", "N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_no1, &
1569 66 : 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1570 : ! 3. TMP_OO2=TMP_OO1.siginv=(T^(tr)FT).siginv
1571 : CALL dbcsr_multiply("N", "N", 1.0_dp, tmp_oo1, almo_scf_env%matrix_sigma_inv(ispin), &
1572 66 : 0.0_dp, tmp_oo2, filter_eps=almo_scf_env%eps_filter)
1573 : ! 4. TMP_OO1=siginv.TMP_OO2=siginv.(T^(tr)FTsiginv)
1574 : CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_sigma_inv(ispin), tmp_oo2, &
1575 66 : 0.0_dp, tmp_oo1, filter_eps=almo_scf_env%eps_filter)
1576 : ! 5. TMP_NO1=T.TMP_OO1.=T.(siginvT^(tr)FTsiginv)
1577 : CALL dbcsr_multiply("N", "N", 1.0_dp, almo_scf_env%matrix_t(ispin), tmp_oo1, &
1578 66 : 0.0_dp, tmp_no1, filter_eps=almo_scf_env%eps_filter)
1579 : ! 6. TMP_NN1=TMP_NO1.T^(tr)=(TsiginvT^(tr)FTsiginv).T^(tr)=RFR
1580 : CALL dbcsr_multiply("N", "T", 1.0_dp, tmp_no1, almo_scf_env%matrix_t(ispin), &
1581 66 : 0.0_dp, tmp_nn1, filter_eps=almo_scf_env%eps_filter)
1582 66 : CALL matrix_almo_to_qs(tmp_nn1, matrix_w(ispin)%matrix, almo_scf_env%mat_distr_aos)
1583 :
1584 66 : CALL dbcsr_release(tmp_nn1)
1585 66 : CALL dbcsr_release(tmp_no1)
1586 66 : CALL dbcsr_release(tmp_oo1)
1587 132 : CALL dbcsr_release(tmp_oo2)
1588 :
1589 : END DO
1590 :
1591 66 : CALL timestop(handle)
1592 :
1593 66 : END SUBROUTINE calculate_w_matrix_almo
1594 :
1595 : END MODULE almo_scf_qs
1596 :
|
