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