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