(git:0de0cc2)
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 ! **************************************************************************************************
20  almo_scf_env_type
22  USE cell_types, ONLY: cell_type,&
23  pbc
24  USE cp_control_types, ONLY: dft_control_type
29  cp_fm_struct_type
30  USE cp_fm_types, ONLY: cp_fm_create,&
31  cp_fm_release,&
32  cp_fm_type
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
52  USE kinds, ONLY: dp
53  USE message_passing, ONLY: mp_comm_type
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,&
62  USE qs_ks_types, ONLY: qs_ks_did_change,&
63  qs_ks_env_type,&
65  USE qs_mo_types, ONLY: allocate_mo_set,&
67  init_mo_set,&
68  mo_set_type
72  neighbor_list_iterator_p_type,&
74  neighbor_list_set_p_type
76  USE qs_rho_types, ONLY: qs_rho_get,&
77  qs_rho_type
78  USE qs_scf_types, ONLY: qs_scf_env_type,&
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, &
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  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)
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)
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)
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 
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 
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 
1604 END MODULE almo_scf_qs
1605 
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
Interface between ALMO SCF and QS.
Definition: almo_scf_qs.F:14
subroutine, public almo_scf_update_ks_energy(qs_env, energy, energy_singles_corr)
update qs_env total energy
Definition: almo_scf_qs.F:811
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
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
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
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
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
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
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
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
Types for all ALMO-based methods.
integer, parameter, public almo_mat_dim_occ
integer, parameter, public almo_mat_dim_virt_full
integer, parameter, public almo_mat_dim_aobasis
integer, parameter, public almo_mat_dim_virt
integer, parameter, public almo_mat_dim_virt_disc
Define the atomic kind types and their sub types.
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.
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
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
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
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
various routines to log and control the output. The idea is that decisions about where to log should ...
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
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
unit conversion facility
Definition: cp_units.F:30
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
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public almo_constraint_distance
integer, parameter, public do_bondparm_covalent
integer, parameter, public almo_mat_distr_molecular
integer, parameter, public almo_domain_layout_molecular
integer, parameter, public do_bondparm_vdw
integer, parameter, public almo_mat_distr_atomic
integer, parameter, public almo_constraint_block_diagonal
integer, parameter, public almo_constraint_ao_overlap
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Interface to the message passing library MPI.
Define the data structure for the molecule information.
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.
Define the data structure for the particle information.
Perform a QUICKSTEP wavefunction optimization (single point)
Definition: qs_energy.F:14
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.
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.
routines that build the Kohn-Sham matrix (i.e calculate the coulomb and xc parts
Definition: qs_ks_methods.F:22
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)
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
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
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
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
subroutine, public deallocate_mo_set(mo_set)
Deallocate a wavefunction data structure.
Definition: qs_mo_types.F:352
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
Define the neighbor list data types and the corresponding functionality.
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
subroutine, public neighbor_list_iterator_release(iterator_set)
...
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
methods of the rho structure (defined in qs_rho_types)
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...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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
module that contains the definitions of the scf types
Definition: qs_scf_types.F:14
subroutine, public scf_env_create(scf_env)
allocates and initialize an scf_env
Definition: qs_scf_types.F:132