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