1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate localized minimal basis and analyze wavefunctions 8!> \par History 9!> 12.2016 created [JGH] 10!> \author JGH 11! ************************************************************************************************** 12MODULE minbas_wfn_analysis 13 USE atomic_charges, ONLY: print_atomic_charges,& 14 print_bond_orders 15 USE atomic_kind_types, ONLY: atomic_kind_type 16 USE bibliography, ONLY: Lu2004,& 17 cite_reference 18 USE cell_types, ONLY: cell_type 19 USE cp_blacs_env, ONLY: cp_blacs_env_type 20 USE cp_control_types, ONLY: dft_control_type 21 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,& 22 cp_dbcsr_plus_fm_fm_t,& 23 cp_dbcsr_sm_fm_multiply,& 24 dbcsr_allocate_matrix_set,& 25 dbcsr_deallocate_matrix_set 26 USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale 27 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 28 cp_fm_struct_release,& 29 cp_fm_struct_type 30 USE cp_fm_types, ONLY: cp_fm_create,& 31 cp_fm_get_diag,& 32 cp_fm_release,& 33 cp_fm_to_fm,& 34 cp_fm_type 35 USE cp_gemm_interface, ONLY: cp_gemm 36 USE cp_log_handling, ONLY: cp_get_default_logger,& 37 cp_logger_type 38 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 39 cp_print_key_unit_nr 40 USE cp_para_types, ONLY: cp_para_env_type 41 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube 42 USE dbcsr_api, ONLY: & 43 dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_dot, dbcsr_get_block_p, & 44 dbcsr_get_occupation, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, & 45 dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, & 46 dbcsr_p_type, dbcsr_release, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry, & 47 dbcsr_type_symmetric 48 USE input_section_types, ONLY: section_get_ivals,& 49 section_vals_get,& 50 section_vals_get_subs_vals,& 51 section_vals_type,& 52 section_vals_val_get 53 USE iterate_matrix, ONLY: invert_Hotelling 54 USE kinds, ONLY: default_path_length,& 55 dp 56 USE message_passing, ONLY: mp_sum 57 USE minbas_methods, ONLY: minbas_calculation 58 USE molden_utils, ONLY: write_mos_molden 59 USE mulliken, ONLY: compute_bond_order,& 60 mulliken_charges 61 USE particle_list_types, ONLY: particle_list_type 62 USE particle_methods, ONLY: get_particle_set 63 USE particle_types, ONLY: particle_type 64 USE pw_env_types, ONLY: pw_env_get,& 65 pw_env_type 66 USE pw_pool_types, ONLY: pw_pool_create_pw,& 67 pw_pool_give_back_pw,& 68 pw_pool_type 69 USE pw_types, ONLY: COMPLEXDATA1D,& 70 REALDATA3D,& 71 REALSPACE,& 72 RECIPROCALSPACE,& 73 pw_p_type 74 USE qs_collocate_density, ONLY: calculate_wavefunction 75 USE qs_environment_types, ONLY: get_qs_env,& 76 qs_environment_type 77 USE qs_kind_types, ONLY: qs_kind_type 78 USE qs_ks_types, ONLY: get_ks_env,& 79 qs_ks_env_type 80 USE qs_mo_methods, ONLY: make_basis_lowdin 81 USE qs_mo_types, ONLY: allocate_mo_set,& 82 deallocate_mo_set,& 83 get_mo_set,& 84 mo_set_p_type,& 85 mo_set_type,& 86 set_mo_set 87 USE qs_subsys_types, ONLY: qs_subsys_get,& 88 qs_subsys_type 89#include "./base/base_uses.f90" 90 91 IMPLICIT NONE 92 PRIVATE 93 94 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minbas_wfn_analysis' 95 96 PUBLIC :: minbas_analysis 97 98! ************************************************************************************************** 99 100CONTAINS 101 102! ************************************************************************************************** 103!> \brief ... 104!> \param qs_env ... 105!> \param input_section ... 106!> \param unit_nr ... 107! ************************************************************************************************** 108 SUBROUTINE minbas_analysis(qs_env, input_section, unit_nr) 109 TYPE(qs_environment_type), POINTER :: qs_env 110 TYPE(section_vals_type), POINTER :: input_section 111 INTEGER, INTENT(IN) :: unit_nr 112 113 CHARACTER(len=*), PARAMETER :: routineN = 'minbas_analysis', & 114 routineP = moduleN//':'//routineN 115 116 INTEGER :: handle, homo, i, ispin, nao, natom, & 117 nimages, nmao, nmo, nspin 118 INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: ecount 119 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, row_blk_sizes 120 LOGICAL :: do_bondorder, explicit, full_ortho, occeq 121 REAL(KIND=dp) :: alpha, amax, eps_filter, filter_eps, & 122 trace 123 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: border, fnorm, mcharge, prmao 124 REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers 125 TYPE(cp_blacs_env_type), POINTER :: blacs_env 126 TYPE(cp_fm_struct_type), POINTER :: fm_struct_a, fm_struct_b, fm_struct_c 127 TYPE(cp_fm_type), POINTER :: fm1, fm2, fm3, fm4, fm_mos 128 TYPE(cp_para_env_type), POINTER :: para_env 129 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 130 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, pqmat, quambo, sqmat 131 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s 132 TYPE(dbcsr_type) :: psmat, sinv, smao, smaox, spmat 133 TYPE(dbcsr_type), POINTER :: smat 134 TYPE(dft_control_type), POINTER :: dft_control 135 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mbas, mos 136 TYPE(mo_set_type), POINTER :: mo_set 137 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 138 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 139 TYPE(qs_ks_env_type), POINTER :: ks_env 140 TYPE(section_vals_type), POINTER :: molden_section 141 142 ! only do MINBAS analysis if explicitly requested 143 CALL section_vals_get(input_section, explicit=explicit) 144 IF (.NOT. explicit) RETURN 145 146 ! k-points? 147 CALL get_qs_env(qs_env, dft_control=dft_control) 148 nspin = dft_control%nspins 149 nimages = dft_control%nimages 150 IF (nimages > 1) THEN 151 IF (unit_nr > 0) THEN 152 WRITE (UNIT=unit_nr, FMT="(T2,A)") & 153 "K-Points: Localized Minimal Basis Analysis not available." 154 END IF 155 END IF 156 IF (nimages > 1) RETURN 157 158 CALL timeset(routineN, handle) 159 160 IF (unit_nr > 0) THEN 161 WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!' 162 WRITE (UNIT=unit_nr, FMT="(T26,A)") "LOCALIZED MINIMAL BASIS ANALYSIS" 163 WRITE (UNIT=unit_nr, FMT="(T18,A)") "W.C. Lu et al, J. Chem. Phys. 120, 2629 (2004)" 164 WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!' 165 END IF 166 CALL cite_reference(Lu2004) 167 168 ! input options 169 CALL section_vals_val_get(input_section, "EPS_FILTER", r_val=eps_filter) 170 CALL section_vals_val_get(input_section, "FULL_ORTHOGONALIZATION", l_val=full_ortho) 171 CALL section_vals_val_get(input_section, "BOND_ORDER", l_val=do_bondorder) 172 173 ! generate MAOs and QUAMBOs 174 CALL get_qs_env(qs_env, mos=mos) 175 NULLIFY (quambo, mao_coef) 176 CALL minbas_calculation(qs_env, mos, quambo, mao=mao_coef, iounit=unit_nr, & 177 full_ortho=full_ortho, eps_filter=eps_filter) 178 IF (ASSOCIATED(quambo)) THEN 179 CALL get_mo_set(mo_set=mos(1)%mo_set, nao=nao, nmo=nmo) 180 CALL get_qs_env(qs_env=qs_env, ks_env=ks_env) 181 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom) 182 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist) 183 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom)) 184 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes) 185 CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes) 186 nmao = SUM(col_blk_sizes) 187 188 NULLIFY (pqmat, sqmat) 189 CALL dbcsr_allocate_matrix_set(sqmat, nspin) 190 CALL dbcsr_allocate_matrix_set(pqmat, nspin) 191 DO ispin = 1, nspin 192 ALLOCATE (sqmat(ispin)%matrix) 193 CALL dbcsr_create(matrix=sqmat(ispin)%matrix, & 194 name="SQMAT", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, & 195 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0) 196 ALLOCATE (pqmat(ispin)%matrix) 197 CALL dbcsr_create(matrix=pqmat(ispin)%matrix, & 198 name="PQMAT", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, & 199 row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0) 200 END DO 201 DEALLOCATE (row_blk_sizes, col_blk_sizes) 202 203 ! Start wfn analysis 204 IF (unit_nr > 0) THEN 205 WRITE (unit_nr, '(/,T2,A)') 'Localized Minimal Basis Wavefunction Analysis' 206 END IF 207 208 ! localization of basis 209 DO ispin = 1, nspin 210 amax = dbcsr_get_occupation(quambo(ispin)%matrix) 211 IF (unit_nr > 0) THEN 212 WRITE (unit_nr, '(/,T2,A,I2,T69,F10.3,A2,/)') & 213 'Occupation of Basis Function Representation (Spin) ', ispin, amax*100._dp, ' %' 214 END IF 215 END DO 216 217 CALL get_qs_env(qs_env, matrix_s_kp=matrix_s) 218 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env) 219 CALL cp_fm_struct_create(fm_struct_a, nrow_global=nao, ncol_global=nmao, & 220 para_env=para_env, context=blacs_env) 221 CALL cp_fm_create(fm1, fm_struct_a) 222 CALL cp_fm_struct_create(fm_struct_b, nrow_global=nmao, ncol_global=nmo, & 223 para_env=para_env, context=blacs_env) 224 CALL cp_fm_create(fm2, fm_struct_b) 225 CALL cp_fm_create(fm3, fm_struct_b) 226 CALL cp_fm_struct_create(fm_struct_c, nrow_global=nmo, ncol_global=nmo, & 227 para_env=para_env, context=blacs_env) 228 CALL cp_fm_create(fm4, fm_struct_c) 229 ALLOCATE (fnorm(nmo, nspin), ecount(natom, 3, nspin), prmao(natom, nspin)) 230 ecount = 0 231 prmao = 0.0_dp 232 DO ispin = 1, nspin 233 CALL dbcsr_create(smao, name="S*QM", template=mao_coef(1)%matrix) 234 smat => matrix_s(1, 1)%matrix 235 CALL dbcsr_multiply("N", "N", 1.0_dp, smat, quambo(ispin)%matrix, 0.0_dp, smao) 236 ! calculate atomic extend of basis 237 CALL pm_extend(quambo(ispin)%matrix, smao, ecount(:, :, ispin)) 238 CALL dbcsr_create(sinv, name="QM*S*QM", template=sqmat(ispin)%matrix) 239 CALL dbcsr_multiply("T", "N", 1.0_dp, quambo(ispin)%matrix, smao, 0.0_dp, sqmat(ispin)%matrix) 240 ! atomic MAO projection 241 CALL project_mao(mao_coef(ispin)%matrix, smao, sqmat(ispin)%matrix, prmao(:, ispin)) 242 ! invert overlap 243 CALL invert_Hotelling(sinv, sqmat(ispin)%matrix, 1.e-6_dp, silent=.TRUE.) 244 CALL dbcsr_create(smaox, name="S*QM*SINV", template=smao) 245 CALL dbcsr_multiply("N", "N", 1.0_dp, smao, sinv, 0.0_dp, smaox) 246 CALL copy_dbcsr_to_fm(smaox, fm1) 247 CALL get_mo_set(mos(ispin)%mo_set, mo_coeff=fm_mos, homo=homo) 248 CALL cp_gemm("T", "N", nmao, nmo, nao, 1.0_dp, fm1, fm_mos, 0.0_dp, fm2) 249 CALL cp_dbcsr_sm_fm_multiply(sqmat(ispin)%matrix, fm2, fm3, nmo) 250 CALL cp_gemm("T", "N", nmo, nmo, nmao, 1.0_dp, fm2, fm3, 0.0_dp, fm4) 251 CALL cp_fm_get_diag(fm4, fnorm(1:nmo, ispin)) 252 ! fm2 are the projected MOs (in MAO basis); orthogonalize the occupied subspace 253 CALL make_basis_lowdin(vmatrix=fm2, ncol=homo, matrix_s=sqmat(ispin)%matrix) 254 ! pmat 255 CALL get_mo_set(mos(ispin)%mo_set, occupation_numbers=occupation_numbers, maxocc=alpha) 256 occeq = ALL(occupation_numbers(1:homo) == alpha) 257 CALL dbcsr_copy(pqmat(ispin)%matrix, sqmat(ispin)%matrix) 258 CALL dbcsr_set(pqmat(ispin)%matrix, 0.0_dp) 259 IF (occeq) THEN 260 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pqmat(ispin)%matrix, matrix_v=fm2, & 261 ncol=homo, alpha=alpha, keep_sparsity=.FALSE.) 262 ELSE 263 CALL cp_fm_to_fm(fm2, fm3) 264 CALL cp_fm_column_scale(fm3, occupation_numbers(1:homo)) 265 alpha = 1.0_dp 266 CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=pqmat(ispin)%matrix, matrix_v=fm2, & 267 matrix_g=fm3, ncol=homo, alpha=alpha, keep_sparsity=.TRUE.) 268 END IF 269 270 CALL dbcsr_release(smao) 271 CALL dbcsr_release(smaox) 272 CALL dbcsr_release(sinv) 273 END DO 274 ! Basis extension 275 CALL mp_sum(ecount, para_env%group) 276 IF (unit_nr > 0) THEN 277 IF (nspin == 1) THEN 278 WRITE (unit_nr, '(T2,A,T20,A,T40,A,T60,A)') 'Ref. Atom', ' # > 0.100 ', ' # > 0.010 ', ' # > 0.001 ' 279 DO i = 1, natom 280 WRITE (unit_nr, '(T2,I8,T20,I10,T40,I10,T60,I10)') i, ecount(i, 1:3, 1) 281 END DO 282 ELSE 283 WRITE (unit_nr, '(T2,A,T20,A,T40,A,T60,A)') 'Ref. Atom', ' # > 0.100 ', ' # > 0.010 ', ' # > 0.001 ' 284 DO i = 1, natom 285 WRITE (unit_nr, '(T2,I8,T20,2I6,T40,2I6,T60,2I6)') & 286 i, ecount(i, 1, 1:2), ecount(i, 2, 1:2), ecount(i, 3, 1:2) 287 END DO 288 END IF 289 END IF 290 ! MAO projection 291 CALL mp_sum(prmao, para_env%group) 292 IF (unit_nr > 0) THEN 293 DO ispin = 1, nspin 294 WRITE (unit_nr, '(/,T2,A,I2)') 'Projection on same atom MAO orbitals: Spin ', ispin 295 DO i = 1, natom, 2 296 IF (i < natom) THEN 297 WRITE (unit_nr, '(T2,A,I8,T20,A,F10.6,T42,A,I8,T60,A,F10.6)') & 298 " Atom:", i, "Projection:", prmao(i, ispin), " Atom:", i + 1, "Projection:", prmao(i + 1, ispin) 299 ELSE 300 WRITE (unit_nr, '(T2,A,I8,T20,A,F10.6)') " Atom:", i, "Projection:", prmao(i, ispin) 301 END IF 302 END DO 303 END DO 304 END IF 305 ! MO expansion completness 306 DO ispin = 1, nspin 307 CALL get_mo_set(mos(ispin)%mo_set, homo=homo, nmo=nmo) 308 IF (unit_nr > 0) THEN 309 WRITE (unit_nr, '(/,T2,A,I2)') 'MO expansion in Localized Minimal Basis: Spin ', ispin 310 WRITE (unit_nr, '(T64,A)') 'Occupied Orbitals' 311 WRITE (unit_nr, '(8F10.6)') fnorm(1:homo, ispin) 312 WRITE (unit_nr, '(T65,A)') 'Virtual Orbitals' 313 WRITE (unit_nr, '(8F10.6)') fnorm(homo + 1:nmo, ispin) 314 END IF 315 END DO 316 ! Mulliken population 317 IF (unit_nr > 0) THEN 318 WRITE (unit_nr, '(/,T2,A)') 'Mulliken Population Analysis ' 319 END IF 320 ALLOCATE (mcharge(natom, nspin)) 321 DO ispin = 1, nspin 322 CALL dbcsr_dot(pqmat(ispin)%matrix, sqmat(ispin)%matrix, trace) 323 IF (unit_nr > 0) THEN 324 WRITE (unit_nr, '(T2,A,I2,T66,F15.4)') 'Number of Electrons: Trace(PS) Spin ', ispin, trace 325 END IF 326 CALL mulliken_charges(pqmat(ispin)%matrix, sqmat(ispin)%matrix, para_env, mcharge(:, ispin)) 327 END DO 328 CALL print_atomic_charges(particle_set, qs_kind_set, unit_nr, "Minimal Basis Mulliken Charges", & 329 electronic_charges=mcharge) 330 ! Mayer bond orders 331 IF (do_bondorder) THEN 332 ALLOCATE (border(natom, natom)) 333 border = 0.0_dp 334 CALL dbcsr_create(psmat, name="PS", template=sqmat(1)%matrix, matrix_type=dbcsr_type_no_symmetry) 335 CALL dbcsr_create(spmat, name="SP", template=sqmat(1)%matrix, matrix_type=dbcsr_type_no_symmetry) 336 filter_eps = 1.e-6_dp 337 DO ispin = 1, nspin 338 CALL dbcsr_multiply("N", "N", 1.0_dp, pqmat(ispin)%matrix, sqmat(ispin)%matrix, 0.0_dp, psmat, & 339 filter_eps=filter_eps) 340 CALL dbcsr_multiply("N", "N", 1.0_dp, sqmat(ispin)%matrix, pqmat(ispin)%matrix, 0.0_dp, spmat, & 341 filter_eps=filter_eps) 342 CALL compute_bond_order(psmat, spmat, border) 343 END DO 344 CALL mp_sum(border, para_env%group) 345 border = border*REAL(nspin, KIND=dp) 346 CALL dbcsr_release(psmat) 347 CALL dbcsr_release(spmat) 348 CALL print_bond_orders(particle_set, unit_nr, border) 349 DEALLOCATE (border) 350 END IF 351 352 ! for printing purposes we now copy the QUAMBOs into MO format 353 ALLOCATE (mbas(nspin)) 354 DO ispin = 1, nspin 355 NULLIFY (mbas(ispin)%mo_set) 356 CALL allocate_mo_set(mbas(ispin)%mo_set, nao, nmao, nmao, 0.0_dp, 1.0_dp, 0.0_dp) 357 mo_set => mbas(ispin)%mo_set 358 CALL set_mo_set(mo_set, homo=nmao) 359 ALLOCATE (mo_set%eigenvalues(nmao)) 360 mo_set%eigenvalues = 0.0_dp 361 ALLOCATE (mo_set%occupation_numbers(nmao)) 362 mo_set%occupation_numbers = 1.0_dp 363 CALL cp_fm_create(mo_set%mo_coeff, fm_struct_a) 364 CALL copy_dbcsr_to_fm(quambo(ispin)%matrix, mo_set%mo_coeff) 365 END DO 366 367 ! Print basis functions: cube files 368 DO ispin = 1, nspin 369 mo_set => mbas(ispin)%mo_set 370 CALL get_mo_set(mo_set, mo_coeff=fm_mos) 371 CALL post_minbas_cubes(qs_env, input_section, fm_mos, ispin) 372 END DO 373 ! Print basis functions: molden format 374 molden_section => section_vals_get_subs_vals(input_section, "MINBAS_MOLDEN") 375 CALL write_mos_molden(mbas, qs_kind_set, particle_set, molden_section) 376 DO ispin = 1, nspin 377 CALL deallocate_mo_set(mbas(ispin)%mo_set) 378 END DO 379 DEALLOCATE (mbas) 380 381 DEALLOCATE (fnorm, ecount, prmao, mcharge) 382 CALL cp_fm_release(fm1) 383 CALL cp_fm_release(fm2) 384 CALL cp_fm_release(fm3) 385 CALL cp_fm_release(fm4) 386 CALL cp_fm_struct_release(fm_struct_a) 387 CALL cp_fm_struct_release(fm_struct_b) 388 CALL cp_fm_struct_release(fm_struct_c) 389 390 ! clean up 391 CALL dbcsr_deallocate_matrix_set(sqmat) 392 CALL dbcsr_deallocate_matrix_set(pqmat) 393 CALL dbcsr_deallocate_matrix_set(mao_coef) 394 CALL dbcsr_deallocate_matrix_set(quambo) 395 396 END IF 397 398 IF (unit_nr > 0) THEN 399 WRITE (unit_nr, '(/,T2,A)') & 400 '!--------------------------END OF MINBAS ANALYSIS-----------------------------!' 401 END IF 402 403 CALL timestop(handle) 404 405 END SUBROUTINE minbas_analysis 406 407! ************************************************************************************************** 408!> \brief ... 409!> \param quambo ... 410!> \param smao ... 411!> \param ecount ... 412! ************************************************************************************************** 413 SUBROUTINE pm_extend(quambo, smao, ecount) 414 TYPE(dbcsr_type) :: quambo, smao 415 INTEGER, DIMENSION(:, :), INTENT(INOUT) :: ecount 416 417 CHARACTER(len=*), PARAMETER :: routineN = 'pm_extend', routineP = moduleN//':'//routineN 418 419 INTEGER :: iatom, jatom, n 420 LOGICAL :: found 421 REAL(KIND=dp) :: wij 422 REAL(KIND=dp), DIMENSION(:, :), POINTER :: qblock, sblock 423 TYPE(dbcsr_iterator_type) :: dbcsr_iter 424 425 CALL dbcsr_iterator_start(dbcsr_iter, quambo) 426 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 427 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, qblock) 428 CALL dbcsr_get_block_p(matrix=smao, row=iatom, col=jatom, BLOCK=sblock, found=found) 429 IF (found) THEN 430 n = SIZE(qblock, 2) 431 wij = ABS(SUM(qblock*sblock))/REAL(n, KIND=dp) 432 IF (wij > 0.1_dp) THEN 433 ecount(jatom, 1) = ecount(jatom, 1) + 1 434 ELSEIF (wij > 0.01_dp) THEN 435 ecount(jatom, 2) = ecount(jatom, 2) + 1 436 ELSEIF (wij > 0.001_dp) THEN 437 ecount(jatom, 3) = ecount(jatom, 3) + 1 438 END IF 439 END IF 440 END DO 441 CALL dbcsr_iterator_stop(dbcsr_iter) 442 443 END SUBROUTINE pm_extend 444 445! ************************************************************************************************** 446!> \brief ... 447!> \param mao ... 448!> \param smao ... 449!> \param sovl ... 450!> \param prmao ... 451! ************************************************************************************************** 452 SUBROUTINE project_mao(mao, smao, sovl, prmao) 453 TYPE(dbcsr_type) :: mao, smao, sovl 454 REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: prmao 455 456 CHARACTER(len=*), PARAMETER :: routineN = 'project_mao', routineP = moduleN//':'//routineN 457 458 INTEGER :: i, iatom, jatom, n 459 LOGICAL :: found 460 REAL(KIND=dp) :: wi 461 REAL(KIND=dp), DIMENSION(:, :), POINTER :: qblock, sblock, so 462 TYPE(dbcsr_iterator_type) :: dbcsr_iter 463 464 CALL dbcsr_iterator_start(dbcsr_iter, mao) 465 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 466 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, qblock) 467 CPASSERT(iatom == jatom) 468 CALL dbcsr_get_block_p(matrix=smao, row=iatom, col=jatom, BLOCK=sblock, found=found) 469 IF (found) THEN 470 CALL dbcsr_get_block_p(matrix=sovl, row=iatom, col=jatom, BLOCK=so, found=found) 471 n = SIZE(qblock, 2) 472 DO i = 1, n 473 wi = SUM(qblock(:, i)*sblock(:, i)) 474 prmao(iatom) = prmao(iatom) + wi/so(i, i) 475 END DO 476 END IF 477 END DO 478 CALL dbcsr_iterator_stop(dbcsr_iter) 479 480 END SUBROUTINE project_mao 481 482! ************************************************************************************************** 483!> \brief Computes and prints the Cube Files for the minimal basis set 484!> \param qs_env ... 485!> \param print_section ... 486!> \param minbas_coeff ... 487!> \param ispin ... 488! ************************************************************************************************** 489 SUBROUTINE post_minbas_cubes(qs_env, print_section, minbas_coeff, ispin) 490 TYPE(qs_environment_type), POINTER :: qs_env 491 TYPE(section_vals_type), POINTER :: print_section 492 TYPE(cp_fm_type), POINTER :: minbas_coeff 493 INTEGER, INTENT(IN) :: ispin 494 495 CHARACTER(len=*), PARAMETER :: routineN = 'post_minbas_cubes', & 496 routineP = moduleN//':'//routineN 497 498 CHARACTER(LEN=default_path_length) :: filename, title 499 INTEGER :: i, i_rep, ivec, iw, j, n_rep, natom 500 INTEGER, DIMENSION(:), POINTER :: blk_sizes, first_bas, ilist, stride 501 LOGICAL :: explicit, mpi_io 502 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 503 TYPE(cell_type), POINTER :: cell 504 TYPE(cp_logger_type), POINTER :: logger 505 TYPE(dft_control_type), POINTER :: dft_control 506 TYPE(particle_list_type), POINTER :: particles 507 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 508 TYPE(pw_env_type), POINTER :: pw_env 509 TYPE(pw_p_type) :: wf_g, wf_r 510 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 511 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 512 TYPE(qs_subsys_type), POINTER :: subsys 513 TYPE(section_vals_type), POINTER :: minbas_section 514 515 minbas_section => section_vals_get_subs_vals(print_section, "MINBAS_CUBE") 516 CALL section_vals_get(minbas_section, explicit=explicit) 517 IF (.NOT. explicit) RETURN 518 519 CPASSERT(ASSOCIATED(minbas_coeff)) 520 521 logger => cp_get_default_logger() 522 stride => section_get_ivals(print_section, "MINBAS_CUBE%STRIDE") 523 524 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, & 525 subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control) 526 CALL qs_subsys_get(subsys, particles=particles) 527 528 CALL get_qs_env(qs_env=qs_env, natom=natom) 529 ALLOCATE (blk_sizes(natom), first_bas(0:natom)) 530 CALL get_particle_set(particle_set, qs_kind_set, nmao=blk_sizes) 531 first_bas(0) = 0 532 DO i = 1, natom 533 first_bas(i) = first_bas(i - 1) + blk_sizes(i) 534 END DO 535 536 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 537 CALL pw_pool_create_pw(auxbas_pw_pool, wf_r%pw, use_data=REALDATA3D, in_space=REALSPACE) 538 CALL pw_pool_create_pw(auxbas_pw_pool, wf_g%pw, use_data=COMPLEXDATA1D, in_space=RECIPROCALSPACE) 539 540 ! loop over list of atoms 541 CALL section_vals_val_get(minbas_section, "ATOM_LIST", n_rep_val=n_rep) 542 IF (n_rep == 0) THEN 543 DO i = 1, natom 544 DO ivec = first_bas(i - 1) + 1, first_bas(i) 545 WRITE (filename, '(a4,I5.5,a1,I1.1)') "MINBAS_", ivec, "_", ispin 546 WRITE (title, *) "MINIMAL BASIS ", ivec, " atom ", i, " spin ", ispin 547 mpi_io = .TRUE. 548 iw = cp_print_key_unit_nr(logger, print_section, "MINBAS_CUBE", extension=".cube", & 549 middle_name=TRIM(filename), file_position="REWIND", log_filename=.FALSE., & 550 mpi_io=mpi_io) 551 CALL calculate_wavefunction(minbas_coeff, ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, & 552 cell, dft_control, particle_set, pw_env) 553 CALL cp_pw_to_cube(wf_r%pw, iw, title, particles=particles, stride=stride, mpi_io=mpi_io) 554 CALL cp_print_key_finished_output(iw, logger, print_section, "MINBAS_CUBE", mpi_io=mpi_io) 555 END DO 556 END DO 557 ELSE 558 DO i_rep = 1, n_rep 559 CALL section_vals_val_get(minbas_section, "ATOM_LIST", i_rep_val=i_rep, i_vals=ilist) 560 DO i = 1, SIZE(ilist, 1) 561 j = ilist(i) 562 DO ivec = first_bas(j - 1) + 1, first_bas(j) 563 WRITE (filename, '(a4,I5.5,a1,I1.1)') "MINBAS_", ivec, "_", ispin 564 WRITE (title, *) "MINIMAL BASIS ", ivec, " atom ", j, " spin ", ispin 565 mpi_io = .TRUE. 566 iw = cp_print_key_unit_nr(logger, print_section, "MINBAS_CUBE", extension=".cube", & 567 middle_name=TRIM(filename), file_position="REWIND", log_filename=.FALSE., & 568 mpi_io=mpi_io) 569 CALL calculate_wavefunction(minbas_coeff, ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, & 570 cell, dft_control, particle_set, pw_env) 571 CALL cp_pw_to_cube(wf_r%pw, iw, title, particles=particles, stride=stride, mpi_io=mpi_io) 572 CALL cp_print_key_finished_output(iw, logger, print_section, "MINBAS_CUBE", mpi_io=mpi_io) 573 END DO 574 END DO 575 END DO 576 END IF 577 DEALLOCATE (blk_sizes, first_bas) 578 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw) 579 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_g%pw) 580 581 END SUBROUTINE post_minbas_cubes 582 583END MODULE minbas_wfn_analysis 584