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