Line data Source code
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 : MODULE qs_fb_env_methods
9 :
10 : USE atomic_kind_types, ONLY: atomic_kind_type,&
11 : get_atomic_kind
12 : USE basis_set_types, ONLY: get_gto_basis_set,&
13 : gto_basis_set_p_type,&
14 : gto_basis_set_type
15 : USE cell_types, ONLY: cell_type
16 : USE cp_blacs_env, ONLY: cp_blacs_env_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: &
19 : dbcsr_create, dbcsr_finalize, dbcsr_get_info, dbcsr_iterator_blocks_left, &
20 : dbcsr_iterator_next_block, dbcsr_iterator_readonly_start, dbcsr_iterator_stop, &
21 : dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_reserve_blocks, &
22 : dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
23 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
24 : dbcsr_allocate_matrix_set,&
25 : dbcsr_deallocate_matrix_set
26 : USE cp_fm_basic_linalg, ONLY: cp_fm_gemm,&
27 : cp_fm_symm,&
28 : cp_fm_triangular_invert,&
29 : cp_fm_triangular_multiply,&
30 : cp_fm_uplo_to_full
31 : USE cp_fm_cholesky, ONLY: cp_fm_cholesky_decompose,&
32 : cp_fm_cholesky_reduce,&
33 : cp_fm_cholesky_restore
34 : USE cp_fm_diag, ONLY: choose_eigv_solver,&
35 : cp_fm_power
36 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
37 : cp_fm_struct_release,&
38 : cp_fm_struct_type
39 : USE cp_fm_types, ONLY: cp_fm_create,&
40 : cp_fm_release,&
41 : cp_fm_set_all,&
42 : cp_fm_to_fm,&
43 : cp_fm_type
44 : USE cp_log_handling, ONLY: cp_get_default_logger,&
45 : cp_logger_type
46 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
47 : cp_print_key_unit_nr
48 : USE cp_units, ONLY: cp_unit_from_cp2k
49 : USE input_constants, ONLY: cholesky_dbcsr,&
50 : cholesky_inverse,&
51 : cholesky_off,&
52 : cholesky_reduce,&
53 : cholesky_restore
54 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
55 : section_vals_type,&
56 : section_vals_val_get
57 : USE kinds, ONLY: default_string_length,&
58 : dp
59 : USE message_passing, ONLY: mp_para_env_type
60 : USE orbital_pointers, ONLY: nco,&
61 : ncoset
62 : USE parallel_gemm_api, ONLY: parallel_gemm
63 : USE particle_types, ONLY: particle_type
64 : USE qs_density_matrices, ONLY: calculate_density_matrix
65 : USE qs_diis, ONLY: qs_diis_b_step
66 : USE qs_environment_types, ONLY: get_qs_env,&
67 : qs_environment_type
68 : USE qs_fb_atomic_halo_types, ONLY: &
69 : fb_atomic_halo_build_halo_atoms, fb_atomic_halo_cost, fb_atomic_halo_create, &
70 : fb_atomic_halo_list_create, fb_atomic_halo_list_nullify, fb_atomic_halo_list_obj, &
71 : fb_atomic_halo_list_set, fb_atomic_halo_list_write_info, &
72 : fb_atomic_halo_nelectrons_estimate_Z, fb_atomic_halo_nullify, fb_atomic_halo_obj, &
73 : fb_atomic_halo_set, fb_atomic_halo_sort, fb_build_pair_radii
74 : USE qs_fb_env_types, ONLY: fb_env_get,&
75 : fb_env_has_data,&
76 : fb_env_obj,&
77 : fb_env_set
78 : USE qs_fb_filter_matrix_methods, ONLY: fb_fltrmat_build,&
79 : fb_fltrmat_build_2
80 : USE qs_fb_trial_fns_types, ONLY: fb_trial_fns_create,&
81 : fb_trial_fns_nullify,&
82 : fb_trial_fns_obj,&
83 : fb_trial_fns_release,&
84 : fb_trial_fns_set
85 : USE qs_integral_utils, ONLY: basis_set_list_setup
86 : USE qs_kind_types, ONLY: get_qs_kind,&
87 : qs_kind_type
88 : USE qs_mo_occupation, ONLY: set_mo_occupation
89 : USE qs_mo_types, ONLY: allocate_mo_set,&
90 : deallocate_mo_set,&
91 : get_mo_set,&
92 : init_mo_set,&
93 : mo_set_type,&
94 : set_mo_set
95 : USE qs_scf_types, ONLY: qs_scf_env_type
96 : USE scf_control_types, ONLY: scf_control_type
97 : USE string_utilities, ONLY: compress,&
98 : uppercase
99 : #include "./base/base_uses.f90"
100 :
101 : IMPLICIT NONE
102 :
103 : PRIVATE
104 :
105 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_fb_env_methods'
106 :
107 : PUBLIC :: fb_env_do_diag, &
108 : fb_env_read_input, &
109 : fb_env_build_rcut_auto, &
110 : fb_env_build_atomic_halos, &
111 : fb_env_write_info
112 :
113 : CONTAINS
114 :
115 : ! **************************************************************************************************
116 : !> \brief Do filtered matrix method diagonalisation
117 : !> \param fb_env : the filter matrix environment
118 : !> \param qs_env : quickstep environment
119 : !> \param matrix_ks : DBCSR system (unfiltered) input KS matrix
120 : !> \param matrix_s : DBCSR system (unfiltered) input overlap matrix
121 : !> \param scf_section : SCF input section
122 : !> \param diis_step : whether we are doing a DIIS step
123 : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
124 : ! **************************************************************************************************
125 80 : SUBROUTINE fb_env_do_diag(fb_env, &
126 : qs_env, &
127 : matrix_ks, &
128 : matrix_s, &
129 : scf_section, &
130 : diis_step)
131 : TYPE(fb_env_obj), INTENT(INOUT) :: fb_env
132 : TYPE(qs_environment_type), POINTER :: qs_env
133 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, matrix_s
134 : TYPE(section_vals_type), POINTER :: scf_section
135 : LOGICAL, INTENT(INOUT) :: diis_step
136 :
137 : CHARACTER(LEN=*), PARAMETER :: routineN = 'fb_env_do_diag'
138 :
139 : CHARACTER(len=2) :: spin_string
140 : CHARACTER(len=default_string_length) :: name
141 : INTEGER :: filtered_nfullrowsORcols_total, handle, homo_filtered, ispin, lfomo_filtered, &
142 : my_nmo, nao, ndep, nelectron, nmo, nmo_filtered, nspin, original_nfullrowsORcols_total
143 80 : INTEGER, DIMENSION(:), POINTER :: filtered_rowORcol_block_sizes, &
144 80 : original_rowORcol_block_sizes
145 : LOGICAL :: collective_com
146 : REAL(kind=dp) :: diis_error, eps_default, eps_diis, eps_eigval, fermi_level, filter_temp, &
147 : flexible_electron_count, KTS_filtered, maxocc, mu_filtered
148 80 : REAL(KIND=dp), DIMENSION(:), POINTER :: eigenvalues, eigenvalues_filtered, occ, &
149 80 : occ_filtered
150 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
151 : TYPE(cp_fm_struct_type), POINTER :: filter_fm_struct, fm_struct
152 : TYPE(cp_fm_type) :: fm_matrix_filter, fm_matrix_filtered_ks, &
153 : fm_matrix_filtered_s, fm_matrix_ortho, &
154 : fm_matrix_work
155 : TYPE(cp_fm_type), POINTER :: mo_coeff, mo_coeff_filtered
156 80 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_filter
157 : TYPE(dbcsr_type) :: matrix_filtered_ks, matrix_filtered_s, &
158 : matrix_tmp
159 : TYPE(dbcsr_type), POINTER :: matrix_filtered_p
160 : TYPE(fb_atomic_halo_list_obj) :: atomic_halos
161 : TYPE(fb_trial_fns_obj) :: trial_fns
162 80 : TYPE(mo_set_type), DIMENSION(:), POINTER :: mos, mos_filtered
163 : TYPE(mp_para_env_type), POINTER :: para_env
164 80 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
165 : TYPE(qs_scf_env_type), POINTER :: scf_env
166 : TYPE(scf_control_type), POINTER :: scf_control
167 :
168 : ! TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: sab_orb
169 :
170 80 : CALL timeset(routineN, handle)
171 :
172 80 : NULLIFY (scf_env, scf_control, para_env, blacs_env, particle_set)
173 80 : NULLIFY (eigenvalues, eigenvalues_filtered, occ, occ_filtered)
174 80 : NULLIFY (mos, mos_filtered)
175 80 : NULLIFY (matrix_filter, matrix_filtered_p)
176 80 : NULLIFY (fm_struct, filter_fm_struct)
177 80 : NULLIFY (mo_coeff_filtered, mo_coeff)
178 : ! NULLIFY(sab_orb)
179 80 : CALL fb_atomic_halo_list_nullify(atomic_halos)
180 80 : CALL fb_trial_fns_nullify(trial_fns)
181 80 : NULLIFY (original_rowORcol_block_sizes, filtered_rowORcol_block_sizes)
182 :
183 : ! get qs_env information
184 : CALL get_qs_env(qs_env=qs_env, &
185 : scf_env=scf_env, &
186 : scf_control=scf_control, &
187 : para_env=para_env, &
188 : blacs_env=blacs_env, &
189 : particle_set=particle_set, &
190 80 : mos=mos)
191 :
192 80 : nspin = SIZE(matrix_ks)
193 :
194 : ! ----------------------------------------------------------------------
195 : ! DIIS step - based on non-filtered matrices and MOs
196 : ! ----------------------------------------------------------------------
197 :
198 160 : DO ispin = 1, nspin
199 : CALL copy_dbcsr_to_fm(matrix_ks(ispin)%matrix, &
200 160 : scf_env%scf_work1(ispin))
201 : END DO
202 :
203 80 : eps_diis = scf_control%eps_diis
204 80 : eps_eigval = EPSILON(0.0_dp)
205 :
206 80 : IF (scf_env%iter_count > 1 .AND. .NOT. scf_env%skip_diis) THEN
207 : CALL qs_diis_b_step(scf_env%scf_diis_buffer, mos, scf_env%scf_work1, &
208 : scf_env%scf_work2, scf_env%iter_delta, &
209 : diis_error, diis_step, eps_diis, scf_control%nmixing, &
210 0 : s_matrix=matrix_s, scf_section=scf_section)
211 : ELSE
212 80 : diis_step = .FALSE.
213 : END IF
214 :
215 80 : IF (diis_step) THEN
216 0 : scf_env%iter_param = diis_error
217 0 : scf_env%iter_method = "DIIS/Filter"
218 : ELSE
219 80 : IF (scf_env%mixing_method == 0) THEN
220 0 : scf_env%iter_method = "NoMix/Filter"
221 80 : ELSE IF (scf_env%mixing_method == 1) THEN
222 0 : scf_env%iter_param = scf_env%p_mix_alpha
223 0 : scf_env%iter_method = "P_Mix/Filter"
224 80 : ELSE IF (scf_env%mixing_method > 1) THEN
225 80 : scf_env%iter_param = scf_env%mixing_store%alpha
226 80 : scf_env%iter_method = TRIM(scf_env%mixing_store%iter_method)//"/Filter"
227 : END IF
228 : END IF
229 :
230 : ! ----------------------------------------------------------------------
231 : ! Construct Filter Matrix
232 : ! ----------------------------------------------------------------------
233 :
234 : CALL fb_env_get(fb_env=fb_env, &
235 : filter_temperature=filter_temp, &
236 : atomic_halos=atomic_halos, &
237 80 : eps_default=eps_default)
238 :
239 : ! construct trial functions
240 80 : CALL get_mo_set(mo_set=mos(1), maxocc=maxocc)
241 80 : CALL fb_env_build_trial_fns_auto(fb_env, qs_env, maxocc)
242 : CALL fb_env_get(fb_env=fb_env, &
243 80 : trial_fns=trial_fns)
244 :
245 : ! allocate filter matrix (matrix_filter(ispin)%matrix are
246 : ! nullified by dbcsr_allocate_matrix_set)
247 80 : CALL dbcsr_allocate_matrix_set(matrix_filter, nspin)
248 160 : DO ispin = 1, nspin
249 : ! get system-wide fermi energy and occupancy, we use this to
250 : ! define the filter function used for the filter matrix
251 : CALL get_mo_set(mo_set=mos(ispin), &
252 : mu=fermi_level, &
253 80 : maxocc=maxocc)
254 : ! get filter matrix name
255 80 : WRITE (spin_string, FMT="(I1)") ispin
256 80 : name = TRIM("FILTER MATRIX SPIN "//spin_string)
257 80 : CALL compress(name)
258 80 : CALL uppercase(name)
259 : ! calculate filter matrix (matrix_s(1) is the overlap, the rest
260 : ! in the array are its derivatives)
261 : CALL fb_env_get(fb_env=fb_env, &
262 80 : collective_com=collective_com)
263 240 : IF (collective_com) THEN
264 : CALL fb_fltrmat_build_2(H_mat=matrix_ks(ispin)%matrix, &
265 : S_mat=matrix_s(1)%matrix, &
266 : atomic_halos=atomic_halos, &
267 : trial_fns=trial_fns, &
268 : para_env=para_env, &
269 : particle_set=particle_set, &
270 : fermi_level=fermi_level, &
271 : filter_temp=filter_temp, &
272 : name=name, &
273 : filter_mat=matrix_filter(ispin)%matrix, &
274 16 : tolerance=eps_default)
275 : ELSE
276 : CALL fb_fltrmat_build(H_mat=matrix_ks(ispin)%matrix, &
277 : S_mat=matrix_s(1)%matrix, &
278 : atomic_halos=atomic_halos, &
279 : trial_fns=trial_fns, &
280 : para_env=para_env, &
281 : particle_set=particle_set, &
282 : fermi_level=fermi_level, &
283 : filter_temp=filter_temp, &
284 : name=name, &
285 : filter_mat=matrix_filter(ispin)%matrix, &
286 64 : tolerance=eps_default)
287 : END IF
288 : END DO ! ispin
289 :
290 : ! ----------------------------------------------------------------------
291 : ! Do Filtered Diagonalisation
292 : ! ----------------------------------------------------------------------
293 :
294 : ! Obtain matrix dimensions. KS and S matrices are symmetric, so
295 : ! row_block_sizes and col_block_sizes should be identical. The
296 : ! same applies to the filtered block sizes. Note that filter
297 : ! matrix will have row_block_sizes equal to that of the original,
298 : ! and col_block_sizes equal to that of the filtered. We assume
299 : ! also that the matrix dimensions are identical for both spin
300 : ! channels.
301 : CALL dbcsr_get_info(matrix_ks(1)%matrix, &
302 : row_blk_size=original_rowORcol_block_sizes, &
303 80 : nfullrows_total=original_nfullrowsORcols_total)
304 : CALL dbcsr_get_info(matrix_filter(1)%matrix, &
305 : col_blk_size=filtered_rowORcol_block_sizes, &
306 80 : nfullcols_total=filtered_nfullrowsORcols_total)
307 :
308 : ! filter diagonalisation works on a smaller basis set, and thus
309 : ! requires a new mo_set (molecular orbitals | eigenvectors) and
310 : ! the corresponding matrix pools for the eigenvector coefficients
311 320 : ALLOCATE (mos_filtered(nspin))
312 160 : DO ispin = 1, nspin
313 : CALL get_mo_set(mo_set=mos(ispin), &
314 : maxocc=maxocc, &
315 : nelectron=nelectron, &
316 80 : flexible_electron_count=flexible_electron_count)
317 : CALL allocate_mo_set(mo_set=mos_filtered(ispin), &
318 : nao=filtered_nfullrowsORcols_total, &
319 : nmo=filtered_nfullrowsORcols_total, &
320 : nelectron=nelectron, &
321 : n_el_f=REAL(nelectron, dp), &
322 : maxocc=maxocc, &
323 160 : flexible_electron_count=flexible_electron_count)
324 : END DO ! ispin
325 :
326 : ! create DBCSR filtered KS matrix, this is reused for each spin
327 : ! channel
328 : ! both row_blk_size and col_blk_size should be that of
329 : ! col_blk_size of the filter matrix
330 : CALL dbcsr_create(matrix=matrix_filtered_ks, template=matrix_ks(1)%matrix, &
331 : name=TRIM("FILTERED_KS_MATRIX"), &
332 : matrix_type=dbcsr_type_no_symmetry, &
333 : row_blk_size=filtered_rowORcol_block_sizes, &
334 80 : col_blk_size=filtered_rowORcol_block_sizes)
335 80 : CALL dbcsr_finalize(matrix_filtered_ks)
336 :
337 : ! create DBCSR filtered S (overlap) matrix. Note that
338 : ! matrix_s(1)%matrix is the original overlap matrix---the rest in
339 : ! the array are derivatives, and it should not depend on
340 : ! spin. HOWEVER, since the filter matrix is constructed from KS
341 : ! matrix, and does depend on spin, the filtered S also becomes
342 : ! spin dependent. Nevertheless this matrix is reused for each spin
343 : ! channel
344 : ! both row_blk_size and col_blk_size should be that of
345 : ! col_blk_size of the filter matrix
346 : CALL dbcsr_create(matrix=matrix_filtered_s, template=matrix_s(1)%matrix, &
347 : name=TRIM("FILTERED_S_MATRIX"), &
348 : matrix_type=dbcsr_type_no_symmetry, &
349 : row_blk_size=filtered_rowORcol_block_sizes, &
350 80 : col_blk_size=filtered_rowORcol_block_sizes)
351 80 : CALL dbcsr_finalize(matrix_filtered_s)
352 :
353 : ! create temporary matrix for constructing filtered KS and S
354 : ! the temporary matrix won't be square
355 : CALL dbcsr_create(matrix=matrix_tmp, template=matrix_s(1)%matrix, &
356 : name=TRIM("TEMPORARY_MATRIX"), &
357 : matrix_type=dbcsr_type_no_symmetry, &
358 : row_blk_size=original_rowORcol_block_sizes, &
359 80 : col_blk_size=filtered_rowORcol_block_sizes)
360 80 : CALL dbcsr_finalize(matrix_tmp)
361 :
362 : ! create fm format matrices used for diagonalisation
363 : CALL cp_fm_struct_create(fmstruct=fm_struct, &
364 : para_env=para_env, &
365 : context=blacs_env, &
366 : nrow_global=filtered_nfullrowsORcols_total, &
367 80 : ncol_global=filtered_nfullrowsORcols_total)
368 : ! both fm_matrix_filtered_s and fm_matrix_filtered_ks are reused
369 : ! for each spin channel
370 : CALL cp_fm_create(fm_matrix_filtered_s, &
371 : fm_struct, &
372 80 : name="FM_MATRIX_FILTERED_S")
373 : CALL cp_fm_create(fm_matrix_filtered_ks, &
374 : fm_struct, &
375 80 : name="FM_MATRIX_FILTERED_KS")
376 : ! creaate work matrix
377 80 : CALL cp_fm_create(fm_matrix_work, fm_struct, name="FM_MATRIX_WORK")
378 80 : CALL cp_fm_create(fm_matrix_ortho, fm_struct, name="FM_MATRIX_ORTHO")
379 : ! all fm matrices are created, so can release fm_struct
380 80 : CALL cp_fm_struct_release(fm_struct)
381 :
382 : ! construct filtered KS, S matrix and diagonalise
383 160 : DO ispin = 1, nspin
384 :
385 : ! construct filtered KS matrix
386 : CALL dbcsr_multiply("N", "N", 1.0_dp, &
387 : matrix_ks(ispin)%matrix, matrix_filter(ispin)%matrix, &
388 80 : 0.0_dp, matrix_tmp)
389 : CALL dbcsr_multiply("T", "N", 1.0_dp, &
390 : matrix_filter(ispin)%matrix, matrix_tmp, &
391 80 : 0.0_dp, matrix_filtered_ks)
392 : ! construct filtered S_matrix
393 : CALL dbcsr_multiply("N", "N", 1.0_dp, &
394 : matrix_s(1)%matrix, matrix_filter(ispin)%matrix, &
395 80 : 0.0_dp, matrix_tmp)
396 : CALL dbcsr_multiply("T", "N", 1.0_dp, &
397 : matrix_filter(ispin)%matrix, matrix_tmp, &
398 80 : 0.0_dp, matrix_filtered_s)
399 :
400 : ! now that we have the filtered KS and S matrices for this spin
401 : ! channel, perform ordinary diagonalisation
402 :
403 : ! convert DBCSR matrices to fm format
404 80 : CALL copy_dbcsr_to_fm(matrix_filtered_s, fm_matrix_filtered_s)
405 80 : CALL copy_dbcsr_to_fm(matrix_filtered_ks, fm_matrix_filtered_ks)
406 :
407 80 : CALL get_mo_set(mos_filtered(ispin), nmo=nmo, nao=nao)
408 :
409 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, &
410 : ncol_global=nmo, para_env=para_env, &
411 80 : context=blacs_env)
412 :
413 : ! setup matrix pools for the molecular orbitals
414 : CALL init_mo_set(mos_filtered(ispin), &
415 : fm_struct=fm_struct, &
416 80 : name="FILTERED_MOS")
417 80 : CALL cp_fm_struct_release(fm_struct)
418 :
419 : ! now diagonalise
420 : CALL fb_env_eigensolver(fm_matrix_filtered_ks, &
421 : fm_matrix_filtered_s, &
422 : mos_filtered(ispin), &
423 : fm_matrix_ortho, &
424 : fm_matrix_work, &
425 : eps_eigval, &
426 : ndep, &
427 240 : scf_env%cholesky_method)
428 : END DO ! ispin
429 :
430 : ! release temporary matrices
431 80 : CALL dbcsr_release(matrix_filtered_s)
432 80 : CALL dbcsr_release(matrix_filtered_ks)
433 80 : CALL cp_fm_release(fm_matrix_filtered_s)
434 80 : CALL cp_fm_release(fm_matrix_filtered_ks)
435 80 : CALL cp_fm_release(fm_matrix_work)
436 80 : CALL cp_fm_release(fm_matrix_ortho)
437 :
438 : ! ----------------------------------------------------------------------
439 : ! Construct New Density Matrix
440 : ! ----------------------------------------------------------------------
441 :
442 : ! calculate filtered molecular orbital occupation numbers and fermi
443 : ! level etc
444 : CALL set_mo_occupation(mo_array=mos_filtered, &
445 80 : smear=scf_control%smear)
446 :
447 : ! get the filtered density matrix and then convert back to the
448 : ! full basis version in scf_env ready to be used outside this
449 : ! subroutine
450 80 : ALLOCATE (matrix_filtered_p)
451 : ! the filtered density matrix should have the same sparse
452 : ! structure as the original density matrix, we must copy the
453 : ! sparse structure here, since construction of the density matrix
454 : ! preserves its sparse form, and therefore matrix_filtered_p must
455 : ! have its blocks allocated here now. We assume the original
456 : ! density matrix scf_env%p_mix_new has the same sparse structure
457 : ! in both spin channels.
458 : CALL dbcsr_create(matrix=matrix_filtered_p, template=scf_env%p_mix_new(1, 1)%matrix, &
459 : name=TRIM("FILTERED_MATRIX_P"), &
460 : row_blk_size=filtered_rowORcol_block_sizes, &
461 80 : col_blk_size=filtered_rowORcol_block_sizes)
462 80 : CALL dbcsr_finalize(matrix_filtered_p)
463 : CALL fb_dbcsr_copy_sparse_struct(matrix_filtered_p, &
464 80 : scf_env%p_mix_new(1, 1)%matrix)
465 : ! old implementation, using sab_orb to allocate the blocks in matrix_filtered_p
466 : ! CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb)
467 : ! CALL cp_dbcsr_alloc_block_from_nbl(matrix_filtered_p, sab_orb)
468 80 : CALL dbcsr_set(matrix_filtered_p, 0.0_dp)
469 :
470 160 : DO ispin = 1, nspin
471 : ! calculate matrix_filtered_p
472 : CALL calculate_density_matrix(mos_filtered(ispin), &
473 80 : matrix_filtered_p)
474 : ! convert back to full basis p
475 : CALL dbcsr_multiply("N", "N", 1.0_dp, &
476 : matrix_filter(ispin)%matrix, matrix_filtered_p, &
477 80 : 0.0_dp, matrix_tmp)
478 : CALL dbcsr_multiply("N", "T", 1.0_dp, &
479 : matrix_tmp, matrix_filter(ispin)%matrix, &
480 : 0.0_dp, scf_env%p_mix_new(ispin, 1)%matrix, &
481 160 : retain_sparsity=.TRUE.)
482 : ! note that we want to retain the sparse structure of
483 : ! scf_env%p_mix_new
484 : END DO ! ispin
485 :
486 : ! release temporary matrices
487 80 : CALL dbcsr_release(matrix_tmp)
488 80 : CALL dbcsr_release(matrix_filtered_p)
489 80 : DEALLOCATE (matrix_filtered_p)
490 :
491 : ! ----------------------------------------------------------------------
492 : ! Update MOs
493 : ! ----------------------------------------------------------------------
494 :
495 : ! we still need to convert mos_filtered back to the full basis
496 : ! version (mos) for this, we need to update mo_coeff (and/or
497 : ! mo_coeff_b --- the DBCSR version, if used) of mos
498 :
499 : ! note also that mo_eigenvalues cannot be fully updated, given
500 : ! that the eigenvalues are computed in a smaller basis, and thus
501 : ! do not give the full spectron. Printing of molecular states
502 : ! (molecular DOS) at each SCF step is therefore not recommended
503 : ! when using this method. The idea is that if one wants a full
504 : ! molecular DOS, then one should perform a full diagonalisation
505 : ! without the filters once the SCF has been achieved.
506 :
507 : ! NOTE: from reading the source code, it appears that mo_coeff_b
508 : ! is actually never used by default (DOUBLE CHECK?!). Even
509 : ! subroutine eigensolver_dbcsr updates mo_coeff, and not
510 : ! mo_coeff_b.
511 :
512 : ! create FM format filter matrix
513 : CALL cp_fm_struct_create(fmstruct=filter_fm_struct, &
514 : para_env=para_env, &
515 : context=blacs_env, &
516 : nrow_global=original_nfullrowsORcols_total, &
517 80 : ncol_global=filtered_nfullrowsORcols_total)
518 : CALL cp_fm_create(fm_matrix_filter, &
519 : filter_fm_struct, &
520 80 : name="FM_MATRIX_FILTER")
521 80 : CALL cp_fm_struct_release(filter_fm_struct)
522 :
523 160 : DO ispin = 1, nspin
524 : ! now the full basis mo_set should only contain the reduced
525 : ! number of eigenvectors and eigenvalues
526 : CALL get_mo_set(mo_set=mos_filtered(ispin), &
527 : homo=homo_filtered, &
528 : lfomo=lfomo_filtered, &
529 : nmo=nmo_filtered, &
530 : eigenvalues=eigenvalues_filtered, &
531 : occupation_numbers=occ_filtered, &
532 : mo_coeff=mo_coeff_filtered, &
533 : kTS=kTS_filtered, &
534 80 : mu=mu_filtered)
535 : ! first set all the relevant scalars
536 : CALL set_mo_set(mo_set=mos(ispin), &
537 : homo=homo_filtered, &
538 : lfomo=lfomo_filtered, &
539 : kTS=kTS_filtered, &
540 80 : mu=mu_filtered)
541 : ! now set the arrays and fm_matrices
542 : CALL get_mo_set(mo_set=mos(ispin), &
543 : nmo=nmo, &
544 : occupation_numbers=occ, &
545 : eigenvalues=eigenvalues, &
546 80 : mo_coeff=mo_coeff)
547 : ! number of mos in original mo_set may sometimes be less than
548 : ! nmo_filtered, so we must make sure we do not go out of bounds
549 80 : my_nmo = MIN(nmo, nmo_filtered)
550 2640 : eigenvalues(:) = 0.0_dp
551 5200 : eigenvalues(1:my_nmo) = eigenvalues_filtered(1:my_nmo)
552 2640 : occ(:) = 0.0_dp
553 5200 : occ(1:my_nmo) = occ_filtered(1:my_nmo)
554 : ! convert mo_coeff_filtered back to original basis
555 80 : CALL cp_fm_set_all(matrix=mo_coeff, alpha=0.0_dp)
556 80 : CALL copy_dbcsr_to_fm(matrix_filter(ispin)%matrix, fm_matrix_filter)
557 : CALL cp_fm_gemm("N", "N", &
558 : original_nfullrowsORcols_total, &
559 : my_nmo, &
560 : filtered_nfullrowsORcols_total, &
561 : 1.0_dp, fm_matrix_filter, mo_coeff_filtered, &
562 320 : 0.0_dp, mo_coeff)
563 :
564 : END DO ! ispin
565 :
566 : ! release temporary matrices
567 80 : CALL cp_fm_release(fm_matrix_filter)
568 :
569 : ! ----------------------------------------------------------------------
570 : ! Final Clean Up
571 : ! ----------------------------------------------------------------------
572 :
573 160 : DO ispin = 1, nspin
574 160 : CALL deallocate_mo_set(mo_set=mos_filtered(ispin))
575 : END DO
576 80 : DEALLOCATE (mos_filtered)
577 80 : CALL dbcsr_deallocate_matrix_set(matrix_filter)
578 :
579 80 : CALL timestop(handle)
580 :
581 480 : END SUBROUTINE fb_env_do_diag
582 :
583 : ! **************************************************************************************************
584 : !> \brief The main parallel eigensolver engine for filter matrix diagonalisation
585 : !> \param fm_KS : the BLACS distributed Kohn-Sham matrix, input only
586 : !> \param fm_S : the BLACS distributed overlap matrix, input only
587 : !> \param mo_set : upon output contains the molecular orbitals (eigenvectors)
588 : !> and eigenvalues
589 : !> \param fm_ortho : one of the work matrices, on output, the BLACS distributed
590 : !> matrix for orthogalising the eigen problem. E.g. if using
591 : !> Cholesky inversse, then the upper triangle part contains
592 : !> the inverse of Cholesky U; if not using Cholesky, then it
593 : !> contains the S^-1/2.
594 : !> \param fm_work : work matrix used by eigen solver
595 : !> \param eps_eigval : used for quenching the small numbers when computing S^-1/2
596 : !> any values less than eps_eigval is truncated to zero.
597 : !> \param ndep : if the overlap is not positive definite, then ndep > 0,
598 : !> and equals to the number of linear dependent basis functions
599 : !> in the filtered basis set
600 : !> \param method : method for solving generalised eigenvalue problem
601 : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
602 : ! **************************************************************************************************
603 160 : SUBROUTINE fb_env_eigensolver(fm_KS, fm_S, mo_set, fm_ortho, &
604 : fm_work, eps_eigval, ndep, method)
605 : TYPE(cp_fm_type), INTENT(IN) :: fm_KS, fm_S
606 : TYPE(mo_set_type), INTENT(IN) :: mo_set
607 : TYPE(cp_fm_type), INTENT(IN) :: fm_ortho, fm_work
608 : REAL(KIND=dp), INTENT(IN) :: eps_eigval
609 : INTEGER, INTENT(OUT) :: ndep
610 : INTEGER, INTENT(IN) :: method
611 :
612 : CHARACTER(len=*), PARAMETER :: routineN = 'fb_env_eigensolver'
613 :
614 : CHARACTER(len=8) :: ndep_string
615 : INTEGER :: handle, info, my_method, nao, nmo
616 80 : REAL(KIND=dp), DIMENSION(:), POINTER :: mo_eigenvalues
617 : TYPE(cp_fm_type), POINTER :: mo_coeff
618 :
619 80 : CALL timeset(routineN, handle)
620 :
621 : CALL get_mo_set(mo_set=mo_set, &
622 : nao=nao, &
623 : nmo=nmo, &
624 : eigenvalues=mo_eigenvalues, &
625 80 : mo_coeff=mo_coeff)
626 80 : my_method = method
627 80 : ndep = 0
628 :
629 : ! first, obtain orthogonalisation (ortho) matrix
630 80 : IF (my_method .NE. cholesky_off) THEN
631 64 : CALL cp_fm_to_fm(fm_S, fm_ortho)
632 64 : CALL cp_fm_cholesky_decompose(fm_ortho, info_out=info)
633 64 : IF (info .NE. 0) THEN
634 : CALL cp_warn(__LOCATION__, &
635 : "Unable to perform Cholesky decomposition on the overlap "// &
636 : "matrix. The new filtered basis may not be linearly "// &
637 : "independent set. Revert to using inverse square-root "// &
638 : "of the overlap. To avoid this warning, you can try "// &
639 0 : "to use a higher filter termperature.")
640 : my_method = cholesky_off
641 : ELSE
642 0 : SELECT CASE (my_method)
643 : CASE (cholesky_dbcsr)
644 : CALL cp_abort(__LOCATION__, &
645 0 : "filter matrix method with CHOLESKY_DBCSR is not yet implemented")
646 : CASE (cholesky_reduce)
647 16 : CALL cp_fm_cholesky_reduce(fm_KS, fm_ortho)
648 16 : CALL choose_eigv_solver(fm_KS, fm_work, mo_eigenvalues)
649 16 : CALL cp_fm_cholesky_restore(fm_work, nmo, fm_ortho, mo_coeff, "SOLVE")
650 : CASE (cholesky_restore)
651 32 : CALL cp_fm_uplo_to_full(fm_KS, fm_work)
652 : CALL cp_fm_cholesky_restore(fm_KS, nao, fm_ortho, fm_work, "SOLVE", &
653 32 : pos="RIGHT")
654 : CALL cp_fm_cholesky_restore(fm_work, nao, fm_ortho, fm_KS, "SOLVE", &
655 32 : pos="LEFT", transa="T")
656 32 : CALL choose_eigv_solver(fm_KS, fm_work, mo_eigenvalues)
657 32 : CALL cp_fm_cholesky_restore(fm_work, nmo, fm_ortho, mo_coeff, "SOLVE")
658 : CASE (cholesky_inverse)
659 16 : CALL cp_fm_triangular_invert(fm_ortho)
660 16 : CALL cp_fm_uplo_to_full(fm_KS, fm_work)
661 : CALL cp_fm_triangular_multiply(fm_ortho, &
662 : fm_KS, &
663 : side="R", &
664 : transpose_tr=.FALSE., &
665 : invert_tr=.FALSE., &
666 : uplo_tr="U", &
667 : n_rows=nao, &
668 : n_cols=nao, &
669 16 : alpha=1.0_dp)
670 : CALL cp_fm_triangular_multiply(fm_ortho, &
671 : fm_KS, &
672 : side="L", &
673 : transpose_tr=.TRUE., &
674 : invert_tr=.FALSE., &
675 : uplo_tr="U", &
676 : n_rows=nao, &
677 : n_cols=nao, &
678 16 : alpha=1.0_dp)
679 16 : CALL choose_eigv_solver(fm_KS, fm_work, mo_eigenvalues)
680 : CALL cp_fm_triangular_multiply(fm_ortho, &
681 : fm_work, &
682 : side="L", &
683 : transpose_tr=.FALSE., &
684 : invert_tr=.FALSE., &
685 : uplo_tr="U", &
686 : n_rows=nao, &
687 : n_cols=nmo, &
688 16 : alpha=1.0_dp)
689 80 : CALL cp_fm_to_fm(fm_work, mo_coeff, nmo, 1, 1)
690 : END SELECT
691 : END IF
692 : END IF
693 : IF (my_method == cholesky_off) THEN
694 : ! calculating ortho as S^-1/2 using diagonalisation of S, and
695 : ! solve accordingly
696 16 : CALL cp_fm_to_fm(fm_S, fm_ortho)
697 : CALL cp_fm_power(fm_ortho, fm_work, -0.5_dp, &
698 16 : eps_eigval, ndep)
699 16 : IF (ndep > 0) THEN
700 0 : WRITE (ndep_string, FMT="(I8)") ndep
701 : CALL cp_warn(__LOCATION__, &
702 0 : "Number of linearly dependent filtered orbitals: "//ndep_string)
703 : END IF
704 : ! solve eigen equatoin using S^-1/2
705 : CALL cp_fm_symm("L", "U", nao, nao, 1.0_dp, fm_KS, fm_ortho, &
706 16 : 0.0_dp, fm_work)
707 : CALL parallel_gemm("T", "N", nao, nao, nao, 1.0_dp, fm_ortho, &
708 16 : fm_work, 0.0_dp, fm_KS)
709 16 : CALL choose_eigv_solver(fm_KS, fm_work, mo_eigenvalues)
710 : CALL parallel_gemm("N", "N", nao, nmo, nao, 1.0_dp, fm_ortho, &
711 16 : fm_work, 0.0_dp, mo_coeff)
712 : END IF
713 :
714 80 : CALL timestop(handle)
715 :
716 80 : END SUBROUTINE fb_env_eigensolver
717 :
718 : ! **************************************************************************************************
719 : !> \brief Read input sections for filter matrix method
720 : !> \param fb_env : the filter matrix environment
721 : !> \param scf_section : SCF input section
722 : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
723 : ! **************************************************************************************************
724 50 : SUBROUTINE fb_env_read_input(fb_env, scf_section)
725 :
726 : TYPE(fb_env_obj), INTENT(INOUT) :: fb_env
727 : TYPE(section_vals_type), POINTER :: scf_section
728 :
729 : CHARACTER(len=*), PARAMETER :: routineN = 'fb_env_read_input'
730 :
731 : INTEGER :: handle
732 : LOGICAL :: l_val
733 : REAL(KIND=dp) :: r_val
734 : TYPE(section_vals_type), POINTER :: fb_section
735 :
736 10 : CALL timeset(routineN, handle)
737 :
738 10 : NULLIFY (fb_section)
739 : fb_section => section_vals_get_subs_vals(scf_section, &
740 10 : "DIAGONALIZATION%FILTER_MATRIX")
741 : ! filter_temperature
742 : CALL section_vals_val_get(fb_section, "FILTER_TEMPERATURE", &
743 10 : r_val=r_val)
744 : CALL fb_env_set(fb_env=fb_env, &
745 10 : filter_temperature=r_val)
746 : ! auto_cutoff_scale
747 : CALL section_vals_val_get(fb_section, "AUTO_CUTOFF_SCALE", &
748 10 : r_val=r_val)
749 : CALL fb_env_set(fb_env=fb_env, &
750 10 : auto_cutoff_scale=r_val)
751 : ! communication model
752 : CALL section_vals_val_get(fb_section, "COLLECTIVE_COMMUNICATION", &
753 10 : l_val=l_val)
754 : CALL fb_env_set(fb_env=fb_env, &
755 10 : collective_com=l_val)
756 : ! eps_default
757 : CALL section_vals_val_get(fb_section, "EPS_FB", &
758 10 : r_val=r_val)
759 : CALL fb_env_set(fb_env=fb_env, &
760 10 : eps_default=r_val)
761 :
762 10 : CALL timestop(handle)
763 :
764 10 : END SUBROUTINE fb_env_read_input
765 :
766 : ! **************************************************************************************************
767 : !> \brief Automatically generate the cutoff radii of atoms used for
768 : !> constructing the atomic halos, based on basis set cutoff
769 : !> ranges for each kind
770 : !> \param fb_env : the filter matrix environment
771 : !> \param qs_env : quickstep environment
772 : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
773 : ! **************************************************************************************************
774 10 : SUBROUTINE fb_env_build_rcut_auto(fb_env, qs_env)
775 : TYPE(fb_env_obj), INTENT(INOUT) :: fb_env
776 : TYPE(qs_environment_type), POINTER :: qs_env
777 :
778 : CHARACTER(len=*), PARAMETER :: routineN = 'fb_env_build_rcut_auto'
779 :
780 : INTEGER :: handle, ikind, nkinds
781 : REAL(KIND=dp) :: auto_cutoff_scale, kind_radius
782 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: rcut
783 : TYPE(dft_control_type), POINTER :: dft_control
784 10 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
785 : TYPE(gto_basis_set_type), POINTER :: basis_set
786 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
787 :
788 10 : CALL timeset(routineN, handle)
789 :
790 10 : NULLIFY (rcut, qs_kind_set, dft_control)
791 :
792 : CALL get_qs_env(qs_env=qs_env, &
793 : qs_kind_set=qs_kind_set, &
794 10 : dft_control=dft_control)
795 : CALL fb_env_get(fb_env=fb_env, &
796 10 : auto_cutoff_scale=auto_cutoff_scale)
797 :
798 10 : nkinds = SIZE(qs_kind_set)
799 30 : ALLOCATE (rcut(nkinds))
800 :
801 : ! reading from the other parts of the code, it seemed that
802 : ! aux_fit_basis_set is only used when do_admm is TRUE. This can be
803 : ! seen from the calls to generate_qs_task_list subroutine in
804 : ! qs_create_task_list, found in qs_environment_methods.F:
805 : ! basis_type is only set as input parameter for do_admm
806 : ! calculations, and if not set, the task list is generated using
807 : ! the default basis_set="ORB".
808 40 : ALLOCATE (basis_set_list(nkinds))
809 10 : IF (dft_control%do_admm) THEN
810 0 : CALL basis_set_list_setup(basis_set_list, "AUX_FIT", qs_kind_set)
811 : ELSE
812 10 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
813 : END IF
814 :
815 20 : DO ikind = 1, nkinds
816 10 : basis_set => basis_set_list(ikind)%gto_basis_set
817 10 : CALL get_gto_basis_set(gto_basis_set=basis_set, kind_radius=kind_radius)
818 20 : rcut(ikind) = kind_radius*auto_cutoff_scale
819 : END DO
820 :
821 : CALL fb_env_set(fb_env=fb_env, &
822 10 : rcut=rcut)
823 :
824 : ! cleanup
825 10 : DEALLOCATE (basis_set_list)
826 :
827 10 : CALL timestop(handle)
828 :
829 20 : END SUBROUTINE fb_env_build_rcut_auto
830 :
831 : ! **************************************************************************************************
832 : !> \brief Builds an fb_atomic_halo_list object using information
833 : !> from fb_env
834 : !> \param fb_env the fb_env object
835 : !> \param qs_env : quickstep environment (need this to access particle)
836 : !> positions and their kinds as well as which particles
837 : !> are local to this process
838 : !> \param scf_section : SCF input section, for printing output
839 : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
840 : ! **************************************************************************************************
841 10 : SUBROUTINE fb_env_build_atomic_halos(fb_env, qs_env, scf_section)
842 : TYPE(fb_env_obj), INTENT(INOUT) :: fb_env
843 : TYPE(qs_environment_type), POINTER :: qs_env
844 : TYPE(section_vals_type), POINTER :: scf_section
845 :
846 : CHARACTER(len=*), PARAMETER :: routineN = 'fb_env_build_atomic_halos'
847 :
848 : INTEGER :: handle, iatom, ihalo, max_natoms_local, natoms_global, natoms_local, nelectrons, &
849 : nhalo_atoms, nkinds_global, owner_id_in_halo
850 10 : INTEGER, DIMENSION(:), POINTER :: halo_atoms, local_atoms
851 : REAL(KIND=dp) :: cost
852 10 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radii
853 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: rcut
854 : TYPE(cell_type), POINTER :: cell
855 : TYPE(fb_atomic_halo_list_obj) :: atomic_halos
856 10 : TYPE(fb_atomic_halo_obj), DIMENSION(:), POINTER :: halos
857 : TYPE(mp_para_env_type), POINTER :: para_env
858 10 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
859 10 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
860 :
861 10 : CALL timeset(routineN, handle)
862 :
863 10 : CPASSERT(fb_env_has_data(fb_env))
864 :
865 10 : NULLIFY (cell, halos, halo_atoms, rcut, particle_set, para_env, &
866 10 : qs_kind_set, local_atoms)
867 10 : CALL fb_atomic_halo_list_nullify(atomic_halos)
868 :
869 : ! get relevant data from fb_env
870 : CALL fb_env_get(fb_env=fb_env, &
871 : rcut=rcut, &
872 : local_atoms=local_atoms, &
873 10 : nlocal_atoms=natoms_local)
874 :
875 : ! create atomic_halos
876 10 : CALL fb_atomic_halo_list_create(atomic_halos)
877 :
878 : ! get the number of atoms and kinds:
879 : CALL get_qs_env(qs_env=qs_env, &
880 : natom=natoms_global, &
881 : particle_set=particle_set, &
882 : qs_kind_set=qs_kind_set, &
883 : nkind=nkinds_global, &
884 : para_env=para_env, &
885 10 : cell=cell)
886 :
887 : ! get the maximum number of local atoms across the procs.
888 10 : max_natoms_local = natoms_local
889 10 : CALL para_env%max(max_natoms_local)
890 :
891 : ! create the halos, one for each local atom
892 70 : ALLOCATE (halos(natoms_local))
893 50 : DO ihalo = 1, natoms_local
894 40 : CALL fb_atomic_halo_nullify(halos(ihalo))
895 50 : CALL fb_atomic_halo_create(halos(ihalo))
896 : END DO
897 : CALL fb_atomic_halo_list_set(atomic_halos=atomic_halos, &
898 : nhalos=natoms_local, &
899 10 : max_nhalos=max_natoms_local)
900 : ! build halos
901 40 : ALLOCATE (pair_radii(nkinds_global, nkinds_global))
902 10 : CALL fb_build_pair_radii(rcut, nkinds_global, pair_radii)
903 10 : ihalo = 0
904 50 : DO iatom = 1, natoms_local
905 40 : ihalo = ihalo + 1
906 : CALL fb_atomic_halo_build_halo_atoms(local_atoms(iatom), &
907 : particle_set, &
908 : cell, &
909 : pair_radii, &
910 : halo_atoms, &
911 : nhalo_atoms, &
912 40 : owner_id_in_halo)
913 : CALL fb_atomic_halo_set(atomic_halo=halos(ihalo), &
914 : owner_atom=local_atoms(iatom), &
915 : owner_id_in_halo=owner_id_in_halo, &
916 : natoms=nhalo_atoms, &
917 40 : halo_atoms=halo_atoms)
918 : ! prepare halo_atoms for another halo, do not deallocate, as
919 : ! original data is being pointed at by the atomic halo data
920 : ! structure
921 40 : NULLIFY (halo_atoms)
922 : ! calculate the number of electrons in each halo
923 : nelectrons = fb_atomic_halo_nelectrons_estimate_Z(halos(ihalo), &
924 40 : particle_set)
925 : ! calculate cost
926 40 : cost = fb_atomic_halo_cost(halos(ihalo), particle_set, qs_kind_set)
927 : CALL fb_atomic_halo_set(atomic_halo=halos(ihalo), &
928 : nelectrons=nelectrons, &
929 40 : cost=cost)
930 : ! sort atomic halo
931 90 : CALL fb_atomic_halo_sort(halos(ihalo))
932 : END DO ! iatom
933 10 : DEALLOCATE (pair_radii)
934 :
935 : ! finalise
936 : CALL fb_atomic_halo_list_set(atomic_halos=atomic_halos, &
937 10 : halos=halos)
938 : CALL fb_env_set(fb_env=fb_env, &
939 10 : atomic_halos=atomic_halos)
940 :
941 : ! print info
942 : CALL fb_atomic_halo_list_write_info(atomic_halos, &
943 : para_env, &
944 10 : scf_section)
945 :
946 10 : CALL timestop(handle)
947 :
948 20 : END SUBROUTINE fb_env_build_atomic_halos
949 :
950 : ! **************************************************************************************************
951 : !> \brief Automatically construct the trial functiosn used for generating
952 : !> the filter matrix. It tries to use the single zeta subset from
953 : !> the system GTO basis set as the trial functions
954 : !> \param fb_env : the filter matrix environment
955 : !> \param qs_env : quickstep environment
956 : !> \param maxocc : maximum occupancy for an orbital
957 : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
958 : ! **************************************************************************************************
959 80 : SUBROUTINE fb_env_build_trial_fns_auto(fb_env, qs_env, maxocc)
960 :
961 : TYPE(fb_env_obj), INTENT(INOUT) :: fb_env
962 : TYPE(qs_environment_type), POINTER :: qs_env
963 : REAL(KIND=dp), INTENT(IN) :: maxocc
964 :
965 : CHARACTER(len=*), PARAMETER :: routineN = 'fb_env_build_trial_fns_auto'
966 :
967 : INTEGER :: counter, handle, icgf, ico, ikind, iset, &
968 : ishell, itrial, lshell, max_n_trial, &
969 : nkinds, nset, old_lshell
970 80 : INTEGER, DIMENSION(:), POINTER :: lmax, nfunctions, nshell
971 80 : INTEGER, DIMENSION(:, :), POINTER :: functions
972 : REAL(KIND=dp) :: zeff
973 : TYPE(dft_control_type), POINTER :: dft_control
974 : TYPE(fb_trial_fns_obj) :: trial_fns
975 80 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
976 : TYPE(gto_basis_set_type), POINTER :: basis_set
977 80 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
978 :
979 80 : CALL timeset(routineN, handle)
980 :
981 80 : CPASSERT(fb_env_has_data(fb_env))
982 80 : NULLIFY (nfunctions, functions, basis_set, basis_set_list, qs_kind_set, dft_control)
983 80 : CALL fb_trial_fns_nullify(trial_fns)
984 :
985 : ! create a new trial_fn object
986 80 : CALL fb_trial_fns_create(trial_fns)
987 :
988 : CALL get_qs_env(qs_env=qs_env, &
989 : qs_kind_set=qs_kind_set, &
990 80 : dft_control=dft_control)
991 :
992 80 : nkinds = SIZE(qs_kind_set)
993 :
994 : ! reading from the other parts of the code, it seemed that
995 : ! aux_fit_basis_set is only used when do_admm is TRUE. This can be
996 : ! seen from the calls to generate_qs_task_list subroutine in
997 : ! qs_create_task_list, found in qs_environment_methods.F:
998 : ! basis_type is only set as input parameter for do_admm
999 : ! calculations, and if not set, the task list is generated using
1000 : ! the default basis_set="ORB".
1001 320 : ALLOCATE (basis_set_list(nkinds))
1002 80 : IF (dft_control%do_admm) THEN
1003 0 : CALL basis_set_list_setup(basis_set_list, "AUX_FIT", qs_kind_set)
1004 : ELSE
1005 80 : CALL basis_set_list_setup(basis_set_list, "ORB", qs_kind_set)
1006 : END IF
1007 :
1008 240 : ALLOCATE (nfunctions(nkinds))
1009 160 : nfunctions = 0
1010 :
1011 160 : DO ikind = 1, nkinds
1012 : ! "gto = gaussian type orbital"
1013 80 : basis_set => basis_set_list(ikind)%gto_basis_set
1014 : CALL get_gto_basis_set(gto_basis_set=basis_set, &
1015 : nset=nset, &
1016 : lmax=lmax, &
1017 80 : nshell=nshell)
1018 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), &
1019 80 : zeff=zeff)
1020 :
1021 160 : bset1: DO iset = 1, nset
1022 : ! old_lshell = lmax(iset)
1023 80 : old_lshell = -1
1024 320 : DO ishell = 1, nshell(iset)
1025 240 : lshell = basis_set%l(ishell, iset)
1026 240 : counter = 0
1027 : ! loop over orbitals within the same l
1028 640 : DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
1029 400 : counter = counter + 1
1030 : ! only include the first zeta orbitals
1031 640 : IF ((lshell .GT. old_lshell) .AND. (counter .LE. nco(lshell))) THEN
1032 320 : nfunctions(ikind) = nfunctions(ikind) + 1
1033 : END IF
1034 : END DO
1035 : ! we have got enough trial functions when we have enough
1036 : ! basis functions to accommodate the number of electrons,
1037 : ! AND that that we have included all the first zeta
1038 : ! orbitals of an angular momentum quantum number l
1039 240 : IF (((lshell .GT. old_lshell) .OR. (lshell .EQ. lmax(iset))) .AND. &
1040 : (maxocc*REAL(nfunctions(ikind), dp) .GE. zeff)) THEN
1041 : EXIT bset1
1042 : END IF
1043 160 : old_lshell = lshell
1044 : END DO
1045 : END DO bset1
1046 : END DO ! ikind
1047 :
1048 : ! now that we have the number of trial functions get the trial
1049 : ! functions
1050 160 : max_n_trial = MAXVAL(nfunctions)
1051 320 : ALLOCATE (functions(max_n_trial, nkinds))
1052 480 : functions(:, :) = 0
1053 : ! redo the loops to get the trial function indices within the basis set
1054 160 : DO ikind = 1, nkinds
1055 : ! "gto = gaussian type orbital"
1056 80 : basis_set => basis_set_list(ikind)%gto_basis_set
1057 : CALL get_gto_basis_set(gto_basis_set=basis_set, &
1058 : nset=nset, &
1059 : lmax=lmax, &
1060 80 : nshell=nshell)
1061 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), &
1062 80 : zeff=zeff)
1063 80 : icgf = 0
1064 80 : itrial = 0
1065 160 : bset2: DO iset = 1, nset
1066 80 : old_lshell = -1
1067 320 : DO ishell = 1, nshell(iset)
1068 240 : lshell = basis_set%l(ishell, iset)
1069 240 : counter = 0
1070 : ! loop over orbitals within the same l
1071 640 : DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
1072 400 : icgf = icgf + 1
1073 400 : counter = counter + 1
1074 : ! only include the first zeta orbitals
1075 640 : IF ((lshell .GT. old_lshell) .AND. (counter .LE. nco(lshell))) THEN
1076 320 : itrial = itrial + 1
1077 320 : functions(itrial, ikind) = icgf
1078 : END IF
1079 : END DO
1080 : ! we have got enough trial functions when we have more
1081 : ! basis functions than the number of electrons (obtained
1082 : ! from atomic z), AND that that we have included all the
1083 : ! first zeta orbitals of an angular momentum quantum
1084 : ! number l
1085 240 : IF (((lshell .GT. old_lshell) .OR. (lshell .EQ. lmax(iset))) .AND. &
1086 : (maxocc*REAL(itrial, dp) .GE. zeff)) THEN
1087 : EXIT bset2
1088 : END IF
1089 160 : old_lshell = lshell
1090 : END DO
1091 : END DO bset2
1092 : END DO ! ikind
1093 :
1094 : ! set trial_functions
1095 : CALL fb_trial_fns_set(trial_fns=trial_fns, &
1096 : nfunctions=nfunctions, &
1097 80 : functions=functions)
1098 : ! set fb_env
1099 : CALL fb_env_set(fb_env=fb_env, &
1100 80 : trial_fns=trial_fns)
1101 80 : CALL fb_trial_fns_release(trial_fns)
1102 :
1103 : ! cleanup
1104 80 : DEALLOCATE (basis_set_list)
1105 :
1106 80 : CALL timestop(handle)
1107 :
1108 80 : END SUBROUTINE fb_env_build_trial_fns_auto
1109 :
1110 : ! **************************************************************************************************
1111 : !> \brief Copy the sparse structure of a DBCSR matrix to another, this
1112 : !> means the other matrix will have the same number of blocks
1113 : !> and their corresponding logical locations allocated, although
1114 : !> the blocks does not have to be the same size as the original
1115 : !> \param matrix_out : DBCSR matrix whose blocks are to be allocated
1116 : !> \param matrix_in : DBCSR matrix with existing sparse structure that
1117 : !> is to be copied
1118 : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
1119 : ! **************************************************************************************************
1120 80 : SUBROUTINE fb_dbcsr_copy_sparse_struct(matrix_out, matrix_in)
1121 :
1122 : TYPE(dbcsr_type), INTENT(INOUT) :: matrix_out
1123 : TYPE(dbcsr_type), INTENT(IN) :: matrix_in
1124 :
1125 : INTEGER :: iatom, jatom, nblkcols_total, &
1126 : nblkrows_total, nblks
1127 80 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cols, rows
1128 : TYPE(dbcsr_iterator_type) :: iter
1129 :
1130 : CALL dbcsr_get_info(matrix=matrix_in, &
1131 : nblkrows_total=nblkrows_total, &
1132 80 : nblkcols_total=nblkcols_total)
1133 :
1134 80 : nblks = nblkrows_total*nblkcols_total
1135 240 : ALLOCATE (rows(nblks))
1136 160 : ALLOCATE (cols(nblks))
1137 5200 : rows(:) = 0
1138 5200 : cols(:) = 0
1139 80 : nblks = 0
1140 80 : CALL dbcsr_iterator_readonly_start(iter, matrix_in)
1141 1520 : DO WHILE (dbcsr_iterator_blocks_left(iter))
1142 1440 : CALL dbcsr_iterator_next_block(iter, iatom, jatom)
1143 1440 : nblks = nblks + 1
1144 1440 : rows(nblks) = iatom
1145 1440 : cols(nblks) = jatom
1146 : END DO
1147 80 : CALL dbcsr_iterator_stop(iter)
1148 80 : CALL dbcsr_reserve_blocks(matrix_out, rows(1:nblks), cols(1:nblks))
1149 80 : CALL dbcsr_finalize(matrix_out)
1150 :
1151 : ! cleanup
1152 80 : DEALLOCATE (rows)
1153 80 : DEALLOCATE (cols)
1154 :
1155 160 : END SUBROUTINE fb_dbcsr_copy_sparse_struct
1156 :
1157 : ! **************************************************************************************************
1158 : !> \brief Write out parameters used for the filter matrix method to
1159 : !> output
1160 : !> \param fb_env : the filter matrix environment
1161 : !> \param qs_env : quickstep environment
1162 : !> \param scf_section : SCF input section
1163 : !> \author Lianheng Tong (LT) lianheng.tong@kcl.ac.uk
1164 : ! **************************************************************************************************
1165 20 : SUBROUTINE fb_env_write_info(fb_env, qs_env, scf_section)
1166 : TYPE(fb_env_obj), INTENT(IN) :: fb_env
1167 : TYPE(qs_environment_type), POINTER :: qs_env
1168 : TYPE(section_vals_type), POINTER :: scf_section
1169 :
1170 : CHARACTER(len=*), PARAMETER :: routineN = 'fb_env_write_info'
1171 :
1172 : CHARACTER(LEN=2) :: element_symbol
1173 : INTEGER :: handle, ikind, nkinds, unit_nr
1174 : LOGICAL :: collective_com
1175 : REAL(KIND=dp) :: auto_cutoff_scale, filter_temperature
1176 10 : REAL(KIND=dp), DIMENSION(:), POINTER :: rcut
1177 10 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1178 : TYPE(cp_logger_type), POINTER :: logger
1179 :
1180 10 : CALL timeset(routineN, handle)
1181 :
1182 10 : NULLIFY (rcut, atomic_kind_set, logger)
1183 :
1184 : CALL get_qs_env(qs_env=qs_env, &
1185 10 : atomic_kind_set=atomic_kind_set)
1186 : CALL fb_env_get(fb_env=fb_env, &
1187 : filter_temperature=filter_temperature, &
1188 : auto_cutoff_scale=auto_cutoff_scale, &
1189 : rcut=rcut, &
1190 10 : collective_com=collective_com)
1191 :
1192 10 : nkinds = SIZE(atomic_kind_set)
1193 :
1194 10 : logger => cp_get_default_logger()
1195 : unit_nr = cp_print_key_unit_nr(logger, scf_section, &
1196 : "PRINT%FILTER_MATRIX", &
1197 10 : extension="")
1198 10 : IF (unit_nr > 0) THEN
1199 5 : IF (collective_com) THEN
1200 : WRITE (UNIT=unit_nr, FMT="(/,A,T71,A)") &
1201 1 : " FILTER_MAT_DIAG| MPI communication method:", &
1202 2 : "Collective"
1203 : ELSE
1204 : WRITE (UNIT=unit_nr, FMT="(/,A,T71,A)") &
1205 4 : " FILTER_MAT_DIAG| MPI communication method:", &
1206 8 : "At each step"
1207 : END IF
1208 : WRITE (UNIT=unit_nr, FMT="(A,T71,g10.4)") &
1209 5 : " FILTER_MAT_DIAG| Filter temperature [K]:", &
1210 10 : cp_unit_from_cp2k(filter_temperature, "K")
1211 : WRITE (UNIT=unit_nr, FMT="(A,T71,f10.4)") &
1212 5 : " FILTER_MAT_DIAG| Filter temperature [a.u.]:", &
1213 10 : filter_temperature
1214 : WRITE (UNIT=unit_nr, FMT="(A,T71,f10.4)") &
1215 5 : " FILTER_MAT_DIAG| Auto atomic cutoff radius scale:", &
1216 10 : auto_cutoff_scale
1217 : WRITE (UNIT=unit_nr, FMT="(A)") &
1218 5 : " FILTER_MAT_DIAG| atomic cutoff radii [a.u.]"
1219 10 : DO ikind = 1, nkinds
1220 : CALL get_atomic_kind(atomic_kind=atomic_kind_set(ikind), &
1221 5 : element_symbol=element_symbol)
1222 : WRITE (UNIT=unit_nr, FMT="(A,A,T71,f10.4)") &
1223 10 : " FILTER_MAT_DIAG| ", element_symbol, rcut(ikind)
1224 : END DO ! ikind
1225 : END IF
1226 : CALL cp_print_key_finished_output(unit_nr, logger, scf_section, &
1227 10 : "PRINT%FILTER_MATRIX")
1228 :
1229 10 : CALL timestop(handle)
1230 :
1231 10 : END SUBROUTINE fb_env_write_info
1232 :
1233 : END MODULE qs_fb_env_methods
|
