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