(git:374b731)
Loading...
Searching...
No Matches
almo_scf_qs.F
Go to the documentation of this file.
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! **************************************************************************************************
14MODULE 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
98CONTAINS
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 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 INTEGER, DIMENSION(:), POINTER :: blk_distr, blk_sizes, block_sizes_new, col_distr_new, &
132 col_sizes_new, distr_new_array, row_distr_new, row_sizes_new
133 LOGICAL :: active, one_dim_is_mo, tr
134 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 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 CALL dbcsr_get_info(matrix_qs, distribution=dist_qs)
163
164 natoms = almo_scf_env%natoms
165 nmols = almo_scf_env%nmolecules
166
167 DO dimen = 1, 2 ! 1 - row, 2 - column dimension
168
169 ! distribution pattern is the same for all matrix types (ao, occ, virt)
170 IF (dimen == 1) THEN !rows
171 CALL dbcsr_distribution_get(dist_qs, row_dist=blk_distr)
172 ELSE !columns
173 CALL dbcsr_distribution_get(dist_qs, col_dist=blk_distr)
174 END IF
175
176 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 CALL dbcsr_get_info(matrix_qs, row_blk_size=blk_sizes)
180
181 ! atomic clustering of AOs
182 IF (almo_scf_env%mat_distr_aos == almo_mat_distr_atomic) THEN
183 ALLOCATE (block_sizes_new(natoms), distr_new_array(natoms))
184 block_sizes_new(:) = blk_sizes(:)
185 distr_new_array(:) = blk_distr(:)
186 ! molecular clustering of AOs
187 ELSE IF (almo_scf_env%mat_distr_aos == almo_mat_distr_molecular) THEN
188 ALLOCATE (block_sizes_new(nmols), distr_new_array(nmols))
189 block_sizes_new(:) = 0
190 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 blk_sizes(iatom)
194 END DO
195 DO imol = 1, nmols
196 distr_new_array(imol) = &
197 blk_distr(almo_scf_env%first_atom_of_domain(imol))
198 END DO
199 ELSE
200 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 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 IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN
212 nlength = natoms
213 ALLOCATE (block_sizes_new(nlength))
214 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 ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN
224 nlength = nmols
225 ALLOCATE (block_sizes_new(nlength))
226 IF (size_keys(dimen) == almo_mat_dim_occ) THEN
227 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 WHERE (block_sizes_new == 0) block_sizes_new = 1
231 ELSE IF (size_keys(dimen) == almo_mat_dim_virt_disc) THEN
232 block_sizes_new(:) = almo_scf_env%nvirt_disc_of_domain(:, spin_key)
233 ELSE IF (size_keys(dimen) == almo_mat_dim_virt_full) THEN
234 block_sizes_new(:) = almo_scf_env%nvirt_full_of_domain(:, spin_key)
235 ELSE IF (size_keys(dimen) == almo_mat_dim_virt) THEN
236 block_sizes_new(:) = almo_scf_env%nvirt_of_domain(:, spin_key)
237 END IF
238 ELSE
239 cpabort("Illegal distribution")
240 END IF
241
242 ELSE
243
244 cpabort("Illegal dimension")
245
246 END IF ! end choosing dim size (occ, virt)
247
248 ! distribution for MOs is copied from AOs
249 ALLOCATE (distr_new_array(nlength))
250 ! atomic clustering
251 IF (almo_scf_env%mat_distr_mos == almo_mat_distr_atomic) THEN
252 distr_new_array(:) = blk_distr(:)
253 ! molecular clustering
254 ELSE IF (almo_scf_env%mat_distr_mos == almo_mat_distr_molecular) THEN
255 DO imol = 1, nmols
256 distr_new_array(imol) = &
257 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 IF (dimen == 1) THEN !rows
264 row_sizes_new => block_sizes_new
265 row_distr_new => distr_new_array
266 ELSE !columns
267 col_sizes_new => block_sizes_new
268 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 reuse_arrays=.true.)
276
277 ! Create the matrix
278 CALL dbcsr_create(matrix_new, name_new, &
279 dist_new, symmetry_new, &
280 row_sizes_new, col_sizes_new, reuse_arrays=.true.)
281 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 IF (init_domains) THEN
286
287 CALL dbcsr_distribution_get(dist_new, mynode=mynode)
288 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 nblkrows_tot = dbcsr_nblkrows_total(matrix_new)
351 DO row = 1, nblkrows_tot
352 tr = .false.
353 iblock_row = row
354 iblock_col = row
355 CALL dbcsr_get_stored_coordinates(matrix_new, iblock_row, iblock_col, hold)
356
357 IF (hold .EQ. mynode) THEN
358
359 active = .true.
360
361 one_dim_is_mo = .false.
362 DO dimen = 1, 2 ! 1 - row, 2 - column dimension
363 IF (size_keys(dimen) == almo_mat_dim_occ) one_dim_is_mo = .true.
364 END DO
365 IF (one_dim_is_mo) THEN
366 IF (almo_scf_env%nocc_of_domain(row, spin_key) == 0) active = .false.
367 END IF
368
369 one_dim_is_mo = .false.
370 DO dimen = 1, 2
371 IF (size_keys(dimen) == almo_mat_dim_virt) one_dim_is_mo = .true.
372 END DO
373 IF (one_dim_is_mo) THEN
374 IF (almo_scf_env%nvirt_of_domain(row, spin_key) == 0) active = .false.
375 END IF
376
377 one_dim_is_mo = .false.
378 DO dimen = 1, 2
379 IF (size_keys(dimen) == almo_mat_dim_virt_disc) one_dim_is_mo = .true.
380 END DO
381 IF (one_dim_is_mo) THEN
382 IF (almo_scf_env%nvirt_disc_of_domain(row, spin_key) == 0) active = .false.
383 END IF
384
385 one_dim_is_mo = .false.
386 DO dimen = 1, 2
387 IF (size_keys(dimen) == almo_mat_dim_virt_full) one_dim_is_mo = .true.
388 END DO
389 IF (one_dim_is_mo) THEN
390 IF (almo_scf_env%nvirt_full_of_domain(row, spin_key) == 0) active = .false.
391 END IF
392
393 IF (active) THEN
394 NULLIFY (p_new_block)
395 CALL dbcsr_reserve_block2d(matrix_new, iblock_row, iblock_col, p_new_block)
396 cpassert(ASSOCIATED(p_new_block))
397 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 CALL dbcsr_finalize(matrix_new)
407
408 CALL timestop(handle)
409
410 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 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 CALL timeset(routinen, handle)
434 !RZK-warning if it's not a N(AO)xN(AO) matrix then stop
435
436 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 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 matrix_type=dbcsr_type_no_symmetry)
445 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 CALL dbcsr_set(matrix_almo, 0.0_dp)
455 CALL dbcsr_complete_redistribute(matrix_qs_nosym, matrix_almo, &
456 keep_sparsity=keep_sparsity);
457 CALL dbcsr_release(matrix_qs_nosym)
458
459 CASE DEFAULT
460 cpabort("")
461 END SELECT
462
463 CALL timestop(handle)
464
465 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 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 CALL timeset(routinen, handle)
485 ! RZK-warning if it's not a N(AO)xN(AO) matrix then stop
486
487 SELECT CASE (mat_distr_aos)
488 CASE (almo_mat_distr_atomic)
489 CALL dbcsr_copy_into_existing(matrix_qs, matrix_almo)
490 CASE (almo_mat_distr_molecular)
491 CALL dbcsr_set(matrix_qs, 0.0_dp)
492 CALL dbcsr_complete_redistribute(matrix_almo, matrix_qs, keep_sparsity=.true.)
493 CASE DEFAULT
494 cpabort("")
495 END SELECT
496
497 CALL timestop(handle)
498
499 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 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 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 POINTER :: sab_orb
525 TYPE(qs_ks_env_type), POINTER :: ks_env
526
527 CALL timeset(routinen, handle)
528
529 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 sab_orb=sab_orb)
538
539 nspin = dft_control%nspins
540
541 ! create matrix_ks in the QS env if necessary
542 IF (.NOT. ASSOCIATED(matrix_qs_ks)) THEN
543 CALL dbcsr_allocate_matrix_set(matrix_qs_ks, nspin)
544 DO ispin = 1, nspin
545 ALLOCATE (matrix_qs_ks(ispin)%matrix)
546 CALL dbcsr_create(matrix_qs_ks(ispin)%matrix, &
547 template=matrix_qs_s(1)%matrix)
548 CALL cp_dbcsr_alloc_block_from_nbl(matrix_qs_ks(ispin)%matrix, sab_orb)
549 CALL dbcsr_set(matrix_qs_ks(ispin)%matrix, 0.0_dp)
550 END DO
551 CALL set_ks_env(ks_env, matrix_ks=matrix_qs_ks)
552 END IF
553
554 ! copy to ALMO
555 DO ispin = 1, nspin
556 CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, &
557 matrix_ks(ispin), mat_distr_aos, .false.)
558 CALL dbcsr_filter(matrix_ks(ispin), eps_filter)
559 END DO
560
561 CALL timestop(handle)
562
563 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 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 TYPE(mo_set_type), DIMENSION(:), POINTER :: mos
585 TYPE(qs_scf_env_type), POINTER :: scf_env
586
587 CALL timeset(routinen, handle)
588
589 ! create and init scf_env (this is necessary to return MOs to qs)
590 NULLIFY (mos, fm_struct_tmp, scf_env)
591 ALLOCATE (scf_env)
592 CALL scf_env_create(scf_env)
593
594 !CALL qs_scf_env_initialize(qs_env, scf_env)
595 CALL set_qs_env(qs_env, scf_env=scf_env)
596 CALL get_qs_env(qs_env, dft_control=dft_control, mos=mos)
597
598 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 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 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 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 para_env=almo_scf_env%para_env)
618
619 CALL cp_fm_create(mo_fm_copy, fm_struct_tmp, name="t_orthogonal_converted_to_fm")
620 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 CALL init_mo_set(mos(ispin), fm_ref=mo_fm_copy, name='fm_mo')
624
625 CALL cp_fm_release(mo_fm_copy)
626
627 END DO
628
629 CALL timestop(handle)
630
631 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
651 TYPE(qs_rho_type), POINTER :: rho
652
653 CALL timeset(routinen, handle)
654
655 NULLIFY (rho, rho_ao)
656 nspins = SIZE(matrix_p)
657 CALL get_qs_env(qs_env, rho=rho)
658 CALL qs_rho_get(rho, rho_ao=rho_ao)
659
660 ! set the new density matrix
661 DO ispin = 1, nspins
662 CALL matrix_almo_to_qs(matrix_p(ispin), &
663 rho_ao(ispin)%matrix, &
664 mat_distr_aos)
665 END DO
666 CALL qs_rho_update_rho(rho, qs_env=qs_env)
667 CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.true.)
668
669 CALL timestop(handle)
670
671 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 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 CALL timeset(routinen, handle)
703
704 IF (PRESENT(smear)) THEN
705 smearing = smear
706 ELSE
707 smearing = .false.
708 END IF
709
710 IF (PRESENT(kts_sum)) THEN
711 entropic_term = kts_sum
712 ELSE
713 entropic_term = 0.0_dp
714 END IF
715
716 NULLIFY (energy)
717 CALL get_qs_env(qs_env, energy=energy)
718 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 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 IF (smearing) THEN
725 energy%total = energy%total - energy%kTS + entropic_term
726 END IF
727
728 energy_total = energy%total
729
730 CALL timestop(handle)
731
732 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_qs_ks
766
767 CALL timeset(routinen, handle)
768
769 IF (PRESENT(smear)) THEN
770 smearing = smear
771 ELSE
772 smearing = .false.
773 END IF
774
775 IF (PRESENT(kts_sum)) THEN
776 entropic_term = kts_sum
777 ELSE
778 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 kts_sum=entropic_term)
785
786 nspins = SIZE(matrix_ks)
787
788 ! get KS matrix from the QS env and convert to the ALMO format
789 CALL get_qs_env(qs_env, matrix_ks=matrix_qs_ks)
790 DO ispin = 1, nspins
791 CALL matrix_qs_to_almo(matrix_qs_ks(ispin)%matrix, &
792 matrix_ks(ispin), &
793 mat_distr_aos, .false.)
794 CALL dbcsr_filter(matrix_ks(ispin), eps_filter)
795 END DO
796
797 CALL timestop(handle)
798
799 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 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 CALL get_qs_env(qs_env, energy=qs_energy)
818
819 IF (PRESENT(energy_singles_corr)) THEN
820 qs_energy%singles_corr = energy_singles_corr
821 ELSE
822 qs_energy%singles_corr = 0.0_dp
823 END IF
824
825 IF (PRESENT(energy)) THEN
826 qs_energy%total = energy
827 END IF
828
829 qs_energy%total = qs_energy%total + qs_energy%singles_corr
830
831 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: current_number_neighbors, domain_entries_cpu, &
856 domain_map_global, domain_map_local, first_atom_of_molecule, global_list, &
857 last_atom_of_molecule, list_length_cpu, list_offset_cpu, local_list, offset_for_cpu
858 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: domain_grid, domain_neighbor_list, &
859 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 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 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
872 TYPE(dbcsr_type) :: matrix_s_sym
873 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
874 TYPE(mp_comm_type) :: group
875 TYPE(neighbor_list_iterator_p_type), &
876 DIMENSION(:), POINTER :: nl_iterator, nl_iterator2
877 TYPE(neighbor_list_set_p_type), DIMENSION(:), &
878 POINTER :: sab_almo
879 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
880
881 CALL timeset(routinen, handle)
882
883 ! get a useful output_unit
884 logger => cp_get_default_logger()
885 IF (logger%para_env%is_source()) THEN
886 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
887 ELSE
888 unit_nr = -1
889 END IF
890
891 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 sab_almo=sab_almo)
899
900 ! if we are dealing with molecules get info about them
901 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 ALLOCATE (first_atom_of_molecule(almo_scf_env%nmolecules))
904 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 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 matrix_type=dbcsr_type_no_symmetry)
914 CALL dbcsr_get_info(almo_scf_env%matrix_s(1), &
915 matrix_type=sym)
916 IF (sym .EQ. dbcsr_type_no_symmetry) THEN
917 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 matrix_s_sym)
921 END IF
922
923 ALLOCATE (almo_scf_env%quench_t(almo_scf_env%nspins))
924 ALLOCATE (almo_scf_env%domain_map(almo_scf_env%nspins))
925
926 !DO ispin=1,almo_scf_env%nspins
927 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 init_domains=.false.)
938
939 ! initialize distance quencher
940 CALL dbcsr_work_create(almo_scf_env%quench_t(ispin), &
941 work_mutable=.true.)
942
943 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 CALL dbcsr_get_info(almo_scf_env%quench_t(ispin), distribution=dist)
947 CALL dbcsr_distribution_get(dist, numnodes=nnodes, group=groupid, mynode=mynode)
948 CALL group%set_handle(groupid)
949
950 ! create global atom neighbor list from the local lists
951 ! first, calculate number of local pairs
952 local_list_length = 0
953 CALL neighbor_list_iterator_create(nl_iterator, sab_almo)
954 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 iatom=iatom2, jatom=jatom2, inode=inode2, nnode=nnode2)
959 !WRITE(*,*) "GET INFO: ",iatom2, jatom2, inode2, nnode2
960 IF (inode2 == 1) THEN
961 local_list_length = local_list_length + nnode2
962 END IF
963 END DO
964 CALL neighbor_list_iterator_release(nl_iterator)
965
966 ! second, extract the local list to an array
967 ALLOCATE (local_list(2*local_list_length))
968 local_list(:) = 0
969 local_list_length = 0
970 CALL neighbor_list_iterator_create(nl_iterator2, sab_almo)
971 DO WHILE (neighbor_list_iterate(nl_iterator2) == 0)
972 CALL get_iterator_info(nl_iterator2, &
973 iatom=iatom2, jatom=jatom2)
974 local_list(2*local_list_length + 1) = iatom2
975 local_list(2*local_list_length + 2) = jatom2
976 local_list_length = local_list_length + 1
977 END DO ! end loop over pairs of atoms
978 CALL neighbor_list_iterator_release(nl_iterator2)
979
980 ! third, communicate local length to the other nodes
981 ALLOCATE (list_length_cpu(nnodes), list_offset_cpu(nnodes))
982 CALL group%allgather(2*local_list_length, list_length_cpu)
983
984 ! fourth, create a global list
985 list_offset_cpu(1) = 0
986 DO inode = 2, nnodes
987 list_offset_cpu(inode) = list_offset_cpu(inode - 1) + &
988 list_length_cpu(inode - 1)
989 END DO
990 global_list_length = list_offset_cpu(nnodes) + list_length_cpu(nnodes)
991
992 ! fifth, communicate all list data
993 ALLOCATE (global_list(global_list_length))
994 CALL group%allgatherv(local_list, global_list, &
995 list_length_cpu, list_offset_cpu)
996 DEALLOCATE (list_length_cpu, list_offset_cpu)
997 DEALLOCATE (local_list)
998
999 ! calculate maximum number of atoms surrounding the domain
1000 ALLOCATE (current_number_neighbors(almo_scf_env%ndomains))
1001 current_number_neighbors(:) = 0
1002 global_list_length = global_list_length/2
1003 DO ipair = 1, global_list_length
1004 iatom2 = global_list(2*(ipair - 1) + 1)
1005 jatom2 = global_list(2*(ipair - 1) + 2)
1006 idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
1007 jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
1008 ! add to the list
1009 current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
1010 ! add j,i with i,j
1011 IF (idomain2 .NE. jdomain2) THEN
1012 current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
1013 END IF
1014 END DO
1015 max_domain_neighbors = maxval(current_number_neighbors)
1016
1017 ! use the global atom neighbor list to create a global domain neighbor list
1018 ALLOCATE (domain_neighbor_list_excessive(ndomains, max_domain_neighbors))
1019 current_number_neighbors(:) = 1
1020 DO ipair = 1, ndomains
1021 domain_neighbor_list_excessive(ipair, 1) = ipair
1022 END DO
1023 DO ipair = 1, global_list_length
1024 iatom2 = global_list(2*(ipair - 1) + 1)
1025 jatom2 = global_list(2*(ipair - 1) + 2)
1026 idomain2 = almo_scf_env%domain_index_of_atom(iatom2)
1027 jdomain2 = almo_scf_env%domain_index_of_atom(jatom2)
1028 already_listed = .false.
1029 DO ineighbor = 1, current_number_neighbors(idomain2)
1030 IF (domain_neighbor_list_excessive(idomain2, ineighbor) .EQ. jdomain2) THEN
1031 already_listed = .true.
1032 EXIT
1033 END IF
1034 END DO
1035 IF (.NOT. already_listed) THEN
1036 ! add to the list
1037 current_number_neighbors(idomain2) = current_number_neighbors(idomain2) + 1
1038 domain_neighbor_list_excessive(idomain2, current_number_neighbors(idomain2)) = jdomain2
1039 ! add j,i with i,j
1040 IF (idomain2 .NE. jdomain2) THEN
1041 current_number_neighbors(jdomain2) = current_number_neighbors(jdomain2) + 1
1042 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 DEALLOCATE (global_list)
1047
1048 max_domain_neighbors = maxval(current_number_neighbors)
1049 ALLOCATE (domain_neighbor_list(ndomains, max_domain_neighbors))
1050 domain_neighbor_list(:, :) = 0
1051 domain_neighbor_list(:, :) = domain_neighbor_list_excessive(:, 1:max_domain_neighbors)
1052 DEALLOCATE (domain_neighbor_list_excessive)
1053
1054 ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1055 ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(max_domain_neighbors*ndomains, 2))
1056 almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1057 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 DO row = 1, nblkrows_tot
1065 DO col = 1, current_number_neighbors(row)
1066 tr = .false.
1067 iblock_row = row
1068 iblock_col = domain_neighbor_list(row, col)
1069 CALL dbcsr_get_stored_coordinates(almo_scf_env%quench_t(ispin), &
1070 iblock_row, iblock_col, hold)
1071
1072 IF (hold .EQ. mynode) THEN
1073
1074 ! Translate indices of distribution blocks to indices of domain blocks
1075 ! Rows are AOs
1076 domain_row = almo_scf_env%domain_index_of_ao_block(iblock_row)
1077 ! Columns are electrons (i.e. MOs)
1078 domain_col = almo_scf_env%domain_index_of_mo_block(iblock_col)
1079
1080 SELECT CASE (almo_scf_env%constraint_type)
1081 CASE (almo_constraint_block_diagonal)
1082
1083 block_active = .false.
1084 ! type of electron groups
1085 IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1086
1087 ! type of ao domains
1088 IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1089
1090 ! ao domains are molecular / electron groups are molecular
1091 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 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 IF (almo_scf_env%domain_layout_aos == almo_domain_layout_molecular) THEN
1106
1107 ! ao domains are molecular / electron groups are atomic
1108 cpabort("Illegal: molecular domains and atomic groups")
1109
1110 ELSE
1111
1112 ! ao domains are atomic / electron groups are atomic
1113 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 IF (block_active) THEN
1122
1123 NULLIFY (p_new_block)
1124 CALL dbcsr_reserve_block2d(almo_scf_env%quench_t(ispin), &
1125 iblock_row, iblock_col, p_new_block)
1126 cpassert(ASSOCIATED(p_new_block))
1127 p_new_block(:, :) = 1.0_dp
1128
1129 IF (domain_map_local_entries .GE. max_domain_neighbors*almo_scf_env%ndomains) THEN
1130 cpabort("weird... max_domain_neighbors is exceeded")
1131 END IF
1132 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1133 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1134 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 IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1142
1143 ! type of ao domains
1144 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 iblock_row, iblock_col, p_new_block, found)
1151 IF (found) THEN
1152 ! CPErrorMessage(cp_failure_level,routineP,"S block not found")
1153 ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
1154 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 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 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 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 iblock_row, iblock_col, p_new_block, found)
1182 IF (found) THEN
1183 ! CPErrorMessage(cp_failure_level,routineP,"S block not found")
1184 ! CPPrecondition(.FALSE.,cp_failure_level,routineP,failure)
1185 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 s0 = -log10(abs(almo_scf_env%quencher_s0))
1195 s1 = -log10(abs(almo_scf_env%quencher_s1))
1196 IF (overlap .EQ. 0.0_dp) THEN
1197 overlap = -log10(abs(almo_scf_env%eps_filter)) + 100.0_dp
1198 ELSE
1199 overlap = -log10(overlap)
1200 END IF
1201 IF (s0 .LT. 0.0_dp) THEN
1202 cpabort("S0 is less than zero")
1203 END IF
1204 IF (s1 .LE. 0.0_dp) THEN
1205 cpabort("S1 is less than or equal to zero")
1206 END IF
1207 IF (s0 .GE. s1) THEN
1208 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 IF (overlap .LT. s1) THEN
1213 NULLIFY (p_new_block)
1214 CALL dbcsr_reserve_block2d(almo_scf_env%quench_t(ispin), &
1215 iblock_row, iblock_col, p_new_block)
1216 cpassert(ASSOCIATED(p_new_block))
1217 IF (overlap .LE. s0) THEN
1218 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 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 IF (abs(p_new_block(1, 1)) .GT. abs(almo_scf_env%eps_filter)) THEN
1228 IF (domain_map_local_entries .GE. max_domain_neighbors*almo_scf_env%ndomains) THEN
1229 cpabort("weird... max_domain_neighbors is exceeded")
1230 END IF
1231 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1232 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1233 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 IF (almo_scf_env%domain_layout_mos == almo_domain_layout_molecular) THEN
1242
1243 ! type of ao domains
1244 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 IF (domain_row == domain_col) THEN
1252 distance = 0.0_dp
1253 contact_atom_1 = first_atom_of_molecule(domain_row)
1254 contact_atom_2 = first_atom_of_molecule(domain_col)
1255 ELSE
1256 distance_squared = 1.0e+100_dp
1257 contact_atom_1 = -1
1258 contact_atom_2 = -1
1259 DO iatom = first_atom_of_molecule(domain_row), last_atom_of_molecule(domain_row)
1260 DO jatom = first_atom_of_molecule(domain_col), last_atom_of_molecule(domain_col)
1261 rab(:) = pbc(particle_set(iatom)%r(:), particle_set(jatom)%r(:), cell)
1262 trial_distance_squared = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
1263 IF (trial_distance_squared .LT. distance_squared) THEN
1264 distance_squared = trial_distance_squared
1265 contact_atom_1 = iatom
1266 contact_atom_2 = jatom
1267 END IF
1268 END DO ! jatom
1269 END DO ! iatom
1270 cpassert(contact_atom_1 .GT. 0)
1271 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 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 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 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 rab(:) = pbc(particle_set(domain_row)%r(:), particle_set(domain_col)%r(:), cell)
1296 distance = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
1297 contact_atom_1 = domain_row
1298 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 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 rcov=contact1_radius)
1308 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1309 rcov=contact2_radius)
1310 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 rvdw=contact1_radius)
1313 CALL get_atomic_kind(atomic_kind=particle_set(contact_atom_2)%atomic_kind, &
1314 rvdw=contact2_radius)
1315 ELSE
1316 cpabort("Illegal quencher_radius_type")
1317 END IF
1318 contact1_radius = cp_unit_to_cp2k(contact1_radius, "angstrom")
1319 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 r0 = almo_scf_env%quencher_r0_factor*(contact1_radius + contact2_radius)
1327 !+almo_scf_env%quencher_r0_shift
1328 r1 = almo_scf_env%quencher_r1_factor*(contact1_radius + contact2_radius)
1329 !+almo_scf_env%quencher_r1_shift
1330
1331 IF (r0 .LT. 0.0_dp) THEN
1332 cpabort("R0 is less than zero")
1333 END IF
1334 IF (r1 .LE. 0.0_dp) THEN
1335 cpabort("R1 is less than or equal to zero")
1336 END IF
1337 IF (r0 .GT. r1) THEN
1338 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 IF (distance .LT. r1) THEN
1343 NULLIFY (p_new_block)
1344 CALL dbcsr_reserve_block2d(almo_scf_env%quench_t(ispin), &
1345 iblock_row, iblock_col, p_new_block)
1346 cpassert(ASSOCIATED(p_new_block))
1347 IF (distance .LE. r0) THEN
1348 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 cpabort("")
1355 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 IF (abs(p_new_block(1, 1)) .GT. abs(almo_scf_env%eps_filter)) THEN
1362 IF (domain_map_local_entries .GE. max_domain_neighbors*almo_scf_env%ndomains) THEN
1363 cpabort("weird... max_domain_neighbors is exceeded")
1364 END IF
1365 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 1) = iblock_row
1366 almo_scf_env%domain_map(ispin)%pairs(domain_map_local_entries + 1, 2) = iblock_col
1367 domain_map_local_entries = domain_map_local_entries + 1
1368 END IF
1369
1370 END IF
1371
1372 CASE DEFAULT
1373 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 DEALLOCATE (domain_neighbor_list)
1382 DEALLOCATE (current_number_neighbors)
1383
1384 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 almo_scf_env%eps_filter)
1389
1390 ! check that both domain_map and quench_t have the same number of entries
1391 nblks = dbcsr_get_num_blocks(almo_scf_env%quench_t(ispin))
1392 IF (nblks .NE. domain_map_local_entries) THEN
1393 cpabort("number of blocks is wrong")
1394 END IF
1395
1396 ! first, communicate map sizes on the other nodes
1397 ALLOCATE (domain_entries_cpu(nnodes), offset_for_cpu(nnodes))
1398 CALL group%allgather(2*domain_map_local_entries, domain_entries_cpu)
1399
1400 ! second, create
1401 offset_for_cpu(1) = 0
1402 DO inode = 2, nnodes
1403 offset_for_cpu(inode) = offset_for_cpu(inode - 1) + &
1404 domain_entries_cpu(inode - 1)
1405 END DO
1406 global_entries = offset_for_cpu(nnodes) + domain_entries_cpu(nnodes)
1407
1408 ! communicate all entries
1409 ALLOCATE (domain_map_global(global_entries))
1410 ALLOCATE (domain_map_local(2*domain_map_local_entries))
1411 DO ientry = 1, domain_map_local_entries
1412 domain_map_local(2*(ientry - 1) + 1) = almo_scf_env%domain_map(ispin)%pairs(ientry, 1)
1413 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 domain_entries_cpu, offset_for_cpu)
1417 DEALLOCATE (domain_entries_cpu, offset_for_cpu)
1418 DEALLOCATE (domain_map_local)
1419
1420 DEALLOCATE (almo_scf_env%domain_map(ispin)%index1)
1421 DEALLOCATE (almo_scf_env%domain_map(ispin)%pairs)
1422 ALLOCATE (almo_scf_env%domain_map(ispin)%index1(ndomains))
1423 ALLOCATE (almo_scf_env%domain_map(ispin)%pairs(global_entries/2, 2))
1424 almo_scf_env%domain_map(ispin)%pairs(:, :) = 0
1425 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 max_neig_loop: DO WHILE (max_neig_fails)
1434 ALLOCATE (domain_grid(almo_scf_env%ndomains, 0:max_neig))
1435 domain_grid(:, :) = 0
1436 ! init the number of collected neighbors
1437 domain_grid(:, 0) = 1
1438 ! loop over the records
1439 global_entries = global_entries/2
1440 DO ientry = 1, global_entries
1441 ! get the center
1442 grid1 = domain_map_global(2*ientry)
1443 ! get the neighbor
1444 ineig = domain_map_global(2*(ientry - 1) + 1)
1445 ! check boundaries
1446 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 DEALLOCATE (domain_grid)
1450 max_neig = max_neig*2
1451 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 DO igrid = 1, domain_grid(grid1, 0)
1457 ! compare the current neighbor with that already in the 'book'
1458 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 neig_temp = ineig
1462 ineig = domain_grid(grid1, igrid)
1463 domain_grid(grid1, igrid) = neig_temp
1464 ELSE
1465 IF (domain_grid(grid1, igrid) .EQ. 0) THEN
1466 ! got the empty slot now - insert the record
1467 domain_grid(grid1, igrid) = ineig
1468 ! increase the record counter but do it outside the loop
1469 delayed_increment = .true.
1470 END IF
1471 END IF
1472 END DO
1473 IF (delayed_increment) THEN
1474 domain_grid(grid1, 0) = domain_grid(grid1, 0) + 1
1475 ELSE
1476 ! should not be here - all records must be inserted
1477 cpabort("all records must be inserted")
1478 END IF
1479 END DO
1480 max_neig_fails = .false.
1481 END DO max_neig_loop
1482 DEALLOCATE (domain_map_global)
1483
1484 ientry = 1
1485 DO idomain = 1, almo_scf_env%ndomains
1486 DO ineig = 1, domain_grid(idomain, 0) - 1
1487 almo_scf_env%domain_map(ispin)%pairs(ientry, 1) = domain_grid(idomain, ineig)
1488 almo_scf_env%domain_map(ispin)%pairs(ientry, 2) = idomain
1489 ientry = ientry + 1
1490 END DO
1491 almo_scf_env%domain_map(ispin)%index1(idomain) = ientry
1492 END DO
1493 DEALLOCATE (domain_grid)
1494
1495 !ENDDO ! ispin
1496 IF (almo_scf_env%nspins .EQ. 2) THEN
1497 CALL dbcsr_copy(almo_scf_env%quench_t(2), &
1498 almo_scf_env%quench_t(1))
1499 almo_scf_env%domain_map(2)%pairs(:, :) = &
1500 almo_scf_env%domain_map(1)%pairs(:, :)
1501 almo_scf_env%domain_map(2)%index1(:) = &
1502 almo_scf_env%domain_map(1)%index1(:)
1503 END IF
1504
1505 CALL dbcsr_release(matrix_s_sym)
1506
1507 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 DEALLOCATE (first_atom_of_molecule)
1510 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 CALL timestop(handle)
1531
1532 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 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 CALL timeset(routinen, handle)
1554
1555 IF (almo_scf_env%nspins == 1) THEN
1556 scaling = 2.0_dp
1557 ELSE
1558 scaling = 1.0_dp
1559 END IF
1560
1561 DO ispin = 1, almo_scf_env%nspins
1562
1563 CALL dbcsr_create(tmp_nn1, template=almo_scf_env%matrix_s(1), &
1564 matrix_type=dbcsr_type_no_symmetry)
1565 CALL dbcsr_create(tmp_no1, template=almo_scf_env%matrix_t(ispin), &
1566 matrix_type=dbcsr_type_no_symmetry)
1567 CALL dbcsr_create(tmp_oo1, template=almo_scf_env%matrix_sigma_inv(ispin), &
1568 matrix_type=dbcsr_type_no_symmetry)
1569 CALL dbcsr_create(tmp_oo2, template=almo_scf_env%matrix_sigma_inv(ispin), &
1570 matrix_type=dbcsr_type_no_symmetry)
1571
1572 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 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 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 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 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 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 0.0_dp, tmp_nn1, filter_eps=almo_scf_env%eps_filter)
1591 CALL matrix_almo_to_qs(tmp_nn1, matrix_w(ispin)%matrix, almo_scf_env%mat_distr_aos)
1592
1593 CALL dbcsr_release(tmp_nn1)
1594 CALL dbcsr_release(tmp_no1)
1595 CALL dbcsr_release(tmp_oo1)
1596 CALL dbcsr_release(tmp_oo2)
1597
1598 END DO
1599
1600 CALL timestop(handle)
1601
1602 END SUBROUTINE calculate_w_matrix_almo
1603
1604END MODULE almo_scf_qs
1605
cell_types::pbc
Definition cell_types.F:92
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition cp_dbcsr_operations.F:81
cp_fm_types::cp_fm_release
Definition cp_fm_types.F:90
almo_scf_qs
Interface between ALMO SCF and QS.
Definition almo_scf_qs.F:14
almo_scf_qs::almo_scf_update_ks_energy
subroutine, public almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr)
update qs_env total energy
Definition almo_scf_qs.F:811
almo_scf_qs::construct_qs_mos
subroutine, public construct_qs_mos(qs_env, almo_scf_env)
Create MOs in the QS env to be able to return ALMOs to QS.
Definition almo_scf_qs.F:574
almo_scf_qs::matrix_qs_to_almo
subroutine, public matrix_qs_to_almo(matrix_qs, matrix_almo, mat_distr_aos, keep_sparsity)
convert between two types of matrices: QS style to ALMO style
Definition almo_scf_qs.F:423
almo_scf_qs::almo_dm_to_almo_ks
subroutine, public almo_dm_to_almo_ks(qs_env, matrix_p, matrix_ks, energy_total, eps_filter, mat_distr_aos, smear, kts_sum)
uses the ALMO density matrix to compute ALMO KS matrix and the new energy
Definition almo_scf_qs.F:752
almo_scf_qs::calculate_w_matrix_almo
subroutine, public calculate_w_matrix_almo(matrix_w, almo_scf_env)
Compute matrix W (energy-weighted density matrix) that is needed for the evaluation of forces.
Definition almo_scf_qs.F:1544
almo_scf_qs::matrix_almo_create
subroutine, public matrix_almo_create(matrix_new, matrix_qs, almo_scf_env, name_new, size_keys, symmetry_new, spin_key, init_domains)
create the ALMO matrix templates
Definition almo_scf_qs.F:117
almo_scf_qs::init_almo_ks_matrix_via_qs
subroutine, public init_almo_ks_matrix_via_qs(qs_env, matrix_ks, mat_distr_aos, eps_filter)
Initialization of the QS and ALMO KS matrix.
Definition almo_scf_qs.F:512
almo_scf_qs::almo_dm_to_qs_env
subroutine, public almo_dm_to_qs_env(qs_env, matrix_p, mat_distr_aos)
return density matrix to the qs_env
Definition almo_scf_qs.F:643
almo_scf_qs::almo_scf_construct_quencher
subroutine, public almo_scf_construct_quencher(qs_env, almo_scf_env)
Creates the matrix that imposes absolute locality on MOs.
Definition almo_scf_qs.F:842
almo_scf_types
Types for all ALMO-based methods.
Definition almo_scf_types.F:15
almo_scf_types::almo_mat_dim_occ
integer, parameter, public almo_mat_dim_occ
Definition almo_scf_types.F:36
almo_scf_types::almo_mat_dim_virt_full
integer, parameter, public almo_mat_dim_virt_full
Definition almo_scf_types.F:36
almo_scf_types::almo_mat_dim_aobasis
integer, parameter, public almo_mat_dim_aobasis
Definition almo_scf_types.F:36
almo_scf_types::almo_mat_dim_virt
integer, parameter, public almo_mat_dim_virt
Definition almo_scf_types.F:36
almo_scf_types::almo_mat_dim_virt_disc
integer, parameter, public almo_mat_dim_virt_disc
Definition almo_scf_types.F:36
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition cp_dbcsr_operations.F:18
cp_fm_struct
represent the structure of a full matrix
Definition cp_fm_struct.F:14
cp_fm_struct::cp_fm_struct_create
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition cp_fm_struct.F:125
cp_fm_struct::cp_fm_struct_release
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition cp_fm_struct.F:320
cp_fm_types
represent a full matrix distributed on many processors
Definition cp_fm_types.F:15
cp_fm_types::cp_fm_create
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition cp_fm_types.F:167
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:567
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_to_cp2k
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1150
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::almo_constraint_distance
integer, parameter, public almo_constraint_distance
Definition input_constants.F:983
input_constants::do_bondparm_covalent
integer, parameter, public do_bondparm_covalent
Definition input_constants.F:412
input_constants::almo_mat_distr_molecular
integer, parameter, public almo_mat_distr_molecular
Definition input_constants.F:980
input_constants::almo_domain_layout_molecular
integer, parameter, public almo_domain_layout_molecular
Definition input_constants.F:976
input_constants::do_bondparm_vdw
integer, parameter, public do_bondparm_vdw
Definition input_constants.F:412
input_constants::almo_mat_distr_atomic
integer, parameter, public almo_mat_distr_atomic
Definition input_constants.F:980
input_constants::almo_constraint_block_diagonal
integer, parameter, public almo_constraint_block_diagonal
Definition input_constants.F:983
input_constants::almo_constraint_ao_overlap
integer, parameter, public almo_constraint_ao_overlap
Definition input_constants.F:983
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule_set_info
subroutine, public get_molecule_set_info(molecule_set, atom_to_mol, mol_to_first_atom, mol_to_last_atom, mol_to_nelectrons, mol_to_nbasis, mol_to_charge, mol_to_multiplicity)
returns information about molecules in the set.
Definition molecule_types.F:454
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
qs_energy_types
Definition qs_energy_types.F:14
qs_energy
Perform a QUICKSTEP wavefunction optimization (single point)
Definition qs_energy.F:14
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::set_qs_env
subroutine, public set_qs_env(qs_env, super_cell, mos, qmmm, qmmm_periodic, ewald_env, ewald_pw, mpools, rho_external, external_vxc, mask, scf_control, rel_control, qs_charges, ks_env, ks_qmmm_env, wf_history, scf_env, active_space, input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, efield, linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, do_transport, transport_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, mp2_env, bs_env, kg_env, force, kpoints, wanniercentres, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Set the QUICKSTEP environment.
Definition qs_environment_types.F:1035
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:502
qs_ks_methods
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition qs_ks_methods.F:22
qs_ks_methods::qs_ks_update_qs_env
subroutine, public qs_ks_update_qs_env(qs_env, calculate_forces, just_energy, print_active)
updates the Kohn Sham matrix of the given qs_env (facility method)
Definition qs_ks_methods.F:1057
qs_ks_types
Definition qs_ks_types.F:15
qs_ks_types::set_ks_env
subroutine, public set_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, kpoints, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env)
...
Definition qs_ks_types.F:567
qs_ks_types::qs_ks_did_change
subroutine, public qs_ks_did_change(ks_env, s_mstruct_changed, rho_changed, potential_changed, full_reset)
tells that some of the things relevant to the ks calculation did change. has to be called when change...
Definition qs_ks_types.F:919
qs_mo_types
Definition and initialisation of the mo data type.
Definition qs_mo_types.F:22
qs_mo_types::allocate_mo_set
subroutine, public allocate_mo_set(mo_set, nao, nmo, nelectron, n_el_f, maxocc, flexible_electron_count)
Allocates a mo set and partially initializes it (nao,nmo,nelectron, and flexible_electron_count are v...
Definition qs_mo_types.F:206
qs_mo_types::deallocate_mo_set
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition qs_mo_types.F:352
qs_mo_types::init_mo_set
subroutine, public init_mo_set(mo_set, fm_pool, fm_ref, fm_struct, name)
initializes an allocated mo_set. eigenvalues, mo_coeff, occupation_numbers are valid only after this ...
Definition qs_mo_types.F:245
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition qs_neighbor_list_types.F:549
qs_rho_methods
methods of the rho structure (defined in qs_rho_types)
Definition qs_rho_methods.F:15
qs_rho_methods::qs_rho_update_rho
subroutine, public qs_rho_update_rho(rho_struct, qs_env, rho_xc_external, local_rho_set, task_list_external, task_list_external_soft, pw_env_external, para_env_external)
updates rho_r and rho_g to the rhorho_ao. if use_kinetic_energy_density also computes tau_r and tau_g...
Definition qs_rho_methods.F:375
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
qs_scf_types
module that contains the definitions of the scf types
Definition qs_scf_types.F:14
qs_scf_types::scf_env_create
subroutine, public scf_env_create(scf_env)
allocates and initialize an scf_env
Definition qs_scf_types.F:132
almo_scf_types::almo_scf_env_type
Definition almo_scf_types.F:104
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_fm_struct::cp_fm_struct_type
keeps the information about the structure of a full matrix
Definition cp_fm_struct.F:82
cp_fm_types::cp_fm_type
represent a full matrix
Definition cp_fm_types.F:116
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_comm_type
Definition message_passing.F:189
molecule_types::molecule_type
Definition molecule_types.F:110
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_ks_types::qs_ks_env_type
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...
Definition qs_ks_types.F:135
qs_mo_types::mo_set_type
Definition qs_mo_types.F:46
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62
qs_scf_types::qs_scf_env_type
Definition qs_scf_types.F:95