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