1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Calculate MAO's and analyze wavefunctions 8!> \par History 9!> 03.2016 created [JGH] 10!> 12.2016 split into four modules [JGH] 11!> \author JGH 12! ************************************************************************************************** 13MODULE mao_wfn_analysis 14 USE atomic_kind_types, ONLY: get_atomic_kind 15 USE basis_set_types, ONLY: gto_basis_set_p_type 16 USE bibliography, ONLY: Ehrhardt1985,& 17 Heinzmann1976,& 18 cite_reference 19 USE cp_blacs_env, ONLY: cp_blacs_env_type 20 USE cp_control_types, ONLY: dft_control_type 21 USE cp_dbcsr_cholesky, ONLY: cp_dbcsr_cholesky_decompose,& 22 cp_dbcsr_cholesky_restore 23 USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl 24 USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,& 25 dbcsr_deallocate_matrix_set 26 USE cp_para_types, ONLY: cp_para_env_type 27 USE dbcsr_api, ONLY: & 28 dbcsr_copy, dbcsr_create, dbcsr_desymmetrize, dbcsr_distribution_type, dbcsr_dot, & 29 dbcsr_get_block_diag, dbcsr_get_block_p, dbcsr_get_info, dbcsr_iterator_blocks_left, & 30 dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, & 31 dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_replicate_all, & 32 dbcsr_reserve_diag_blocks, dbcsr_type, dbcsr_type_no_symmetry, dbcsr_type_symmetric 33 USE input_section_types, ONLY: section_vals_get,& 34 section_vals_type,& 35 section_vals_val_get 36 USE iterate_matrix, ONLY: invert_Hotelling 37 USE kinds, ONLY: dp 38 USE kpoint_types, ONLY: kpoint_type 39 USE mao_methods, ONLY: mao_basis_analysis,& 40 mao_build_q,& 41 mao_reference_basis 42 USE mao_optimizer, ONLY: mao_optimize 43 USE mathlib, ONLY: invmat_symm 44 USE message_passing, ONLY: mp_sum 45 USE particle_methods, ONLY: get_particle_set 46 USE particle_types, ONLY: particle_type 47 USE qs_environment_types, ONLY: get_qs_env,& 48 qs_environment_type 49 USE qs_kind_types, ONLY: get_qs_kind,& 50 qs_kind_type 51 USE qs_ks_types, ONLY: get_ks_env,& 52 qs_ks_env_type 53 USE qs_neighbor_list_types, ONLY: get_iterator_info,& 54 neighbor_list_iterate,& 55 neighbor_list_iterator_create,& 56 neighbor_list_iterator_p_type,& 57 neighbor_list_iterator_release,& 58 neighbor_list_set_p_type,& 59 release_neighbor_list_sets 60 USE qs_neighbor_lists, ONLY: setup_neighbor_list 61 USE qs_overlap, ONLY: build_overlap_matrix_simple 62 USE qs_rho_types, ONLY: qs_rho_get,& 63 qs_rho_type 64#include "./base/base_uses.f90" 65 66 IMPLICIT NONE 67 PRIVATE 68 69 TYPE block_type 70 REAL(KIND=dp), DIMENSION(:, :), ALLOCATABLE :: mat 71 END TYPE block_type 72 73 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mao_wfn_analysis' 74 75 PUBLIC :: mao_analysis 76 77! ************************************************************************************************** 78 79CONTAINS 80 81! ************************************************************************************************** 82!> \brief ... 83!> \param qs_env ... 84!> \param input_section ... 85!> \param unit_nr ... 86! ************************************************************************************************** 87 SUBROUTINE mao_analysis(qs_env, input_section, unit_nr) 88 TYPE(qs_environment_type), POINTER :: qs_env 89 TYPE(section_vals_type), POINTER :: input_section 90 INTEGER, INTENT(IN) :: unit_nr 91 92 CHARACTER(len=*), PARAMETER :: routineN = 'mao_analysis', routineP = moduleN//':'//routineN 93 94 CHARACTER(len=2) :: element_symbol, esa, esb, esc 95 INTEGER :: fall, handle, ia, iab, iabc, iatom, ib, ic, icol, ikind, irow, ispin, jatom, & 96 mao_basis, max_iter, me, na, nab, nabc, natom, nb, nc, nimages, nspin, ssize 97 INTEGER, DIMENSION(:), POINTER :: col_blk_sizes, mao_blk, mao_blk_sizes, & 98 orb_blk, row_blk_sizes 99 LOGICAL :: analyze_ua, explicit, fo, for, fos, & 100 found, neglect_abc, print_basis 101 REAL(KIND=dp) :: deltaq, electra(2), eps_ab, eps_abc, eps_filter, eps_fun, eps_grad, epsx, & 102 senabc, senmax, threshold, total_charge, total_spin, ua_charge(2), zeff 103 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: occnumA, occnumABC, qab, qmatab, qmatac, & 104 qmatbc, raq, sab, selnABC, sinv, & 105 smatab, smatac, smatbc, uaq 106 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: occnumAB, selnAB 107 REAL(KIND=dp), DIMENSION(:, :), POINTER :: block, cmao, diag, qblka, qblkb, qblkc, & 108 rblkl, rblku, sblk, sblka, sblkb, sblkc 109 TYPE(block_type), ALLOCATABLE, DIMENSION(:) :: rowblock 110 TYPE(cp_blacs_env_type), POINTER :: blacs_env 111 TYPE(cp_para_env_type), POINTER :: para_env 112 TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist 113 TYPE(dbcsr_iterator_type) :: dbcsr_iter 114 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: mao_coef, mao_dmat, mao_qmat, mao_smat, & 115 matrix_q, matrix_smm, matrix_smo 116 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_ks, matrix_p, matrix_s 117 TYPE(dbcsr_type) :: amat, axmat, cgmat, cholmat, crumat, & 118 qmat, qmat_diag, rumat, smat_diag, & 119 sumat, tmat 120 TYPE(dft_control_type), POINTER :: dft_control 121 TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: mao_basis_set_list, orb_basis_set_list 122 TYPE(kpoint_type), POINTER :: kpoints 123 TYPE(neighbor_list_iterator_p_type), & 124 DIMENSION(:), POINTER :: nl_iterator 125 TYPE(neighbor_list_set_p_type), DIMENSION(:), & 126 POINTER :: sab_all, sab_orb, smm_list, smo_list 127 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 128 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 129 TYPE(qs_ks_env_type), POINTER :: ks_env 130 TYPE(qs_rho_type), POINTER :: rho 131 132! only do MAO analysis if explicitely requested 133 134 CALL section_vals_get(input_section, explicit=explicit) 135 IF (.NOT. explicit) RETURN 136 137 CALL timeset(routineN, handle) 138 139 IF (unit_nr > 0) THEN 140 WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!' 141 WRITE (UNIT=unit_nr, FMT="(T36,A)") "MAO ANALYSIS" 142 WRITE (UNIT=unit_nr, FMT="(T12,A)") "Claus Ehrhardt and Reinhart Ahlrichs, TCA 68:231-245 (1985)" 143 WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!' 144 END IF 145 CALL cite_reference(Heinzmann1976) 146 CALL cite_reference(Ehrhardt1985) 147 148 ! input options 149 CALL section_vals_val_get(input_section, "REFERENCE_BASIS", i_val=mao_basis) 150 CALL section_vals_val_get(input_section, "EPS_FILTER", r_val=eps_filter) 151 CALL section_vals_val_get(input_section, "EPS_FUNCTION", r_val=eps_fun) 152 CALL section_vals_val_get(input_section, "EPS_GRAD", r_val=eps_grad) 153 CALL section_vals_val_get(input_section, "MAX_ITER", i_val=max_iter) 154 CALL section_vals_val_get(input_section, "PRINT_BASIS", l_val=print_basis) 155 CALL section_vals_val_get(input_section, "NEGLECT_ABC", l_val=neglect_abc) 156 CALL section_vals_val_get(input_section, "AB_THRESHOLD", r_val=eps_ab) 157 CALL section_vals_val_get(input_section, "ABC_THRESHOLD", r_val=eps_abc) 158 CALL section_vals_val_get(input_section, "ANALYZE_UNASSIGNED_CHARGE", l_val=analyze_ua) 159 160 ! k-points? 161 CALL get_qs_env(qs_env, dft_control=dft_control) 162 nimages = dft_control%nimages 163 IF (nimages > 1) THEN 164 IF (unit_nr > 0) THEN 165 WRITE (UNIT=unit_nr, FMT="(T2,A)") & 166 "K-Points: MAO's determined and analyzed using Gamma-Point only." 167 END IF 168 END IF 169 170 ! Reference basis set 171 NULLIFY (mao_basis_set_list, orb_basis_set_list) 172 CALL mao_reference_basis(qs_env, mao_basis, mao_basis_set_list, orb_basis_set_list, & 173 unit_nr, print_basis) 174 175 ! neighbor lists 176 NULLIFY (smm_list, smo_list) 177 CALL setup_neighbor_list(smm_list, mao_basis_set_list, qs_env=qs_env) 178 CALL setup_neighbor_list(smo_list, mao_basis_set_list, orb_basis_set_list, qs_env=qs_env) 179 180 ! overlap matrices 181 NULLIFY (matrix_smm, matrix_smo) 182 CALL get_qs_env(qs_env, ks_env=ks_env) 183 CALL build_overlap_matrix_simple(ks_env, matrix_smm, & 184 mao_basis_set_list, mao_basis_set_list, smm_list) 185 CALL build_overlap_matrix_simple(ks_env, matrix_smo, & 186 mao_basis_set_list, orb_basis_set_list, smo_list) 187 188 ! get reference density matrix and overlap matrix 189 CALL get_qs_env(qs_env, rho=rho, matrix_s_kp=matrix_s) 190 CALL qs_rho_get(rho, rho_ao_kp=matrix_p) 191 nspin = SIZE(matrix_p, 1) 192 ! 193 ! Q matrix 194 IF (nimages == 1) THEN 195 CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter) 196 ELSE 197 CALL get_qs_env(qs_env, matrix_ks_kp=matrix_ks, kpoints=kpoints) 198 CALL mao_build_q(matrix_q, matrix_p, matrix_s, matrix_smm, matrix_smo, smm_list, electra, eps_filter, & 199 nimages=nimages, kpoints=kpoints, matrix_ks=matrix_ks, sab_orb=sab_orb) 200 END IF 201 202 ! check for extended basis sets 203 fall = 0 204 CALL neighbor_list_iterator_create(nl_iterator, smm_list) 205 DO WHILE (neighbor_list_iterate(nl_iterator) == 0) 206 CALL get_iterator_info(nl_iterator, iatom=iatom, jatom=jatom) 207 IF (iatom <= jatom) THEN 208 irow = iatom 209 icol = jatom 210 ELSE 211 irow = jatom 212 icol = iatom 213 END IF 214 CALL dbcsr_get_block_p(matrix=matrix_p(1, 1)%matrix, & 215 row=irow, col=icol, block=block, found=found) 216 IF (.NOT. found) fall = fall + 1 217 END DO 218 CALL neighbor_list_iterator_release(nl_iterator) 219 220 CALL get_qs_env(qs_env=qs_env, para_env=para_env) 221 CALL mp_sum(fall, para_env%group) 222 IF (unit_nr > 0 .AND. fall > 0) THEN 223 WRITE (UNIT=unit_nr, FMT="(/,T2,A,/,T2,A,/)") & 224 "Warning: Extended MAO basis used with original basis filtered density matrix", & 225 "Warning: Possible errors can be controlled with EPS_PGF_ORB" 226 END IF 227 228 ! MAO matrices 229 CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom) 230 CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist) 231 NULLIFY (mao_coef) 232 CALL dbcsr_allocate_matrix_set(mao_coef, nspin) 233 ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom)) 234 CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes, & 235 basis=mao_basis_set_list) 236 CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes) 237 ! check if MAOs have been specified 238 DO iab = 1, natom 239 IF (col_blk_sizes(iab) < 0) & 240 CPABORT("Number of MAOs has to be specified in KIND section for all elements") 241 END DO 242 DO ispin = 1, nspin 243 ! coeficients 244 ALLOCATE (mao_coef(ispin)%matrix) 245 CALL dbcsr_create(matrix=mao_coef(ispin)%matrix, & 246 name="MAO_COEF", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, & 247 row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, nze=0) 248 CALL dbcsr_reserve_diag_blocks(matrix=mao_coef(ispin)%matrix) 249 END DO 250 DEALLOCATE (row_blk_sizes, col_blk_sizes) 251 252 ! optimize MAOs 253 epsx = 1000.0_dp 254 CALL mao_optimize(mao_coef, matrix_q, matrix_smm, electra, max_iter, eps_grad, epsx, unit_nr) 255 256 ! Analyze the MAO basis 257 CALL mao_basis_analysis(mao_coef, matrix_smm, mao_basis_set_list, particle_set, & 258 qs_kind_set, unit_nr, para_env) 259 260 ! Calculate the overlap and density matrix in the new MAO basis 261 NULLIFY (mao_dmat, mao_smat, mao_qmat) 262 CALL dbcsr_allocate_matrix_set(mao_qmat, nspin) 263 CALL dbcsr_allocate_matrix_set(mao_dmat, nspin) 264 CALL dbcsr_allocate_matrix_set(mao_smat, nspin) 265 CALL dbcsr_get_info(mao_coef(1)%matrix, col_blk_size=col_blk_sizes, distribution=dbcsr_dist) 266 DO ispin = 1, nspin 267 ALLOCATE (mao_dmat(ispin)%matrix) 268 CALL dbcsr_create(mao_dmat(ispin)%matrix, name="MAO density", dist=dbcsr_dist, & 269 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, & 270 col_blk_size=col_blk_sizes, nze=0) 271 ALLOCATE (mao_smat(ispin)%matrix) 272 CALL dbcsr_create(mao_smat(ispin)%matrix, name="MAO overlap", dist=dbcsr_dist, & 273 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, & 274 col_blk_size=col_blk_sizes, nze=0) 275 ALLOCATE (mao_qmat(ispin)%matrix) 276 CALL dbcsr_create(mao_qmat(ispin)%matrix, name="MAO covar density", dist=dbcsr_dist, & 277 matrix_type=dbcsr_type_symmetric, row_blk_size=col_blk_sizes, & 278 col_blk_size=col_blk_sizes, nze=0) 279 END DO 280 CALL dbcsr_create(amat, name="MAO overlap", template=mao_dmat(1)%matrix) 281 CALL dbcsr_create(tmat, name="MAO Overlap Inverse", template=amat) 282 CALL dbcsr_create(qmat, name="MAO covar density", template=amat) 283 CALL dbcsr_create(cgmat, name="TEMP matrix", template=mao_coef(1)%matrix) 284 CALL dbcsr_create(axmat, name="TEMP", template=amat, matrix_type=dbcsr_type_no_symmetry) 285 DO ispin = 1, nspin 286 ! calculate MAO overlap matrix 287 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_smm(1)%matrix, mao_coef(ispin)%matrix, & 288 0.0_dp, cgmat) 289 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, 0.0_dp, amat) 290 ! calculate inverse of MAO overlap 291 threshold = 1.e-8_dp 292 CALL invert_Hotelling(tmat, amat, threshold, norm_convergence=1.e-4_dp, silent=.TRUE.) 293 CALL dbcsr_copy(mao_smat(ispin)%matrix, amat) 294 ! calculate q-matrix q = C*Q*C 295 CALL dbcsr_multiply("N", "N", 1.0_dp, matrix_q(ispin)%matrix, mao_coef(ispin)%matrix, & 296 0.0_dp, cgmat, filter_eps=eps_filter) 297 CALL dbcsr_multiply("T", "N", 1.0_dp, mao_coef(ispin)%matrix, cgmat, & 298 0.0_dp, qmat, filter_eps=eps_filter) 299 CALL dbcsr_copy(mao_qmat(ispin)%matrix, qmat) 300 ! calculate density matrix 301 CALL dbcsr_multiply("N", "N", 1.0_dp, qmat, tmat, 0.0_dp, axmat, filter_eps=eps_filter) 302 CALL dbcsr_multiply("N", "N", 1.0_dp, tmat, axmat, 0.0_dp, mao_dmat(ispin)%matrix, & 303 filter_eps=eps_filter) 304 END DO 305 CALL dbcsr_release(amat) 306 CALL dbcsr_release(tmat) 307 CALL dbcsr_release(qmat) 308 CALL dbcsr_release(cgmat) 309 CALL dbcsr_release(axmat) 310 311 ! calculate unassigned charge : n - Tr PS 312 DO ispin = 1, nspin 313 CALL dbcsr_dot(mao_dmat(ispin)%matrix, mao_smat(ispin)%matrix, ua_charge(ispin)) 314 ua_charge(ispin) = electra(ispin) - ua_charge(ispin) 315 END DO 316 IF (unit_nr > 0) THEN 317 WRITE (unit_nr, *) 318 DO ispin = 1, nspin 319 WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") & 320 "Unassigned charge", "Spin ", ispin, "delta charge =", ua_charge(ispin) 321 END DO 322 END IF 323 324 ! occupation numbers: single atoms 325 ! We use S_A = 1 326 ! At the gamma point we use an effective MIC 327 CALL get_qs_env(qs_env, natom=natom) 328 ALLOCATE (occnumA(natom, nspin)) 329 occnumA = 0.0_dp 330 DO ispin = 1, nspin 331 DO iatom = 1, natom 332 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, & 333 row=iatom, col=iatom, block=block, found=found) 334 IF (found) THEN 335 DO iab = 1, SIZE(block, 1) 336 occnumA(iatom, ispin) = occnumA(iatom, ispin) + block(iab, iab) 337 END DO 338 END IF 339 END DO 340 END DO 341 CALL mp_sum(occnumA, para_env%group) 342 343 ! occupation numbers: atom pairs 344 ALLOCATE (occnumAB(natom, natom, nspin)) 345 occnumAB = 0.0_dp 346 DO ispin = 1, nspin 347 CALL dbcsr_create(qmat_diag, name="MAO diagonal density", template=mao_dmat(1)%matrix) 348 CALL dbcsr_create(smat_diag, name="MAO diagonal overlap", template=mao_dmat(1)%matrix) 349 ! replicate the diagonal blocks of the density and overlap matrices 350 CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag) 351 CALL dbcsr_replicate_all(qmat_diag) 352 CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag) 353 CALL dbcsr_replicate_all(smat_diag) 354 DO ia = 1, natom 355 DO ib = ia + 1, natom 356 iab = 0 357 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, & 358 row=ia, col=ib, block=block, found=found) 359 IF (found) iab = 1 360 CALL mp_sum(iab, para_env%group) 361 CPASSERT(iab <= 1) 362 IF (iab == 0 .AND. para_env%ionode) THEN 363 ! AB block is not available N_AB = N_A + N_B 364 ! Do this only on the "source" processor 365 occnumAB(ia, ib, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin) 366 occnumAB(ib, ia, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin) 367 ELSE IF (found) THEN 368 ! owner of AB block performs calculation 369 na = SIZE(block, 1) 370 nb = SIZE(block, 2) 371 nab = na + nb 372 ALLOCATE (sab(nab, nab), qab(nab, nab), sinv(nab, nab)) 373 ! qmat 374 qab(1:na, na + 1:nab) = block(1:na, 1:nb) 375 qab(na + 1:nab, 1:na) = TRANSPOSE(block(1:na, 1:nb)) 376 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=diag, found=fo) 377 CPASSERT(fo) 378 qab(1:na, 1:na) = diag(1:na, 1:na) 379 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=diag, found=fo) 380 CPASSERT(fo) 381 qab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb) 382 ! smat 383 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, & 384 row=ia, col=ib, block=block, found=fo) 385 CPASSERT(fo) 386 sab(1:na, na + 1:nab) = block(1:na, 1:nb) 387 sab(na + 1:nab, 1:na) = TRANSPOSE(block(1:na, 1:nb)) 388 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=diag, found=fo) 389 CPASSERT(fo) 390 sab(1:na, 1:na) = diag(1:na, 1:na) 391 CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=diag, found=fo) 392 CPASSERT(fo) 393 sab(na + 1:nab, na + 1:nab) = diag(1:nb, 1:nb) 394 ! inv smat 395 sinv(1:nab, 1:nab) = sab(1:nab, 1:nab) 396 CALL invmat_symm(sinv) 397 ! Tr(Q*Sinv) 398 occnumAB(ia, ib, ispin) = SUM(qab*sinv) 399 occnumAB(ib, ia, ispin) = occnumAB(ia, ib, ispin) 400 ! 401 DEALLOCATE (sab, qab, sinv) 402 END IF 403 END DO 404 END DO 405 CALL dbcsr_release(qmat_diag) 406 CALL dbcsr_release(smat_diag) 407 END DO 408 CALL mp_sum(occnumAB, para_env%group) 409 410 ! calculate shared electron numbers (AB) 411 ALLOCATE (selnAB(natom, natom, nspin)) 412 selnAB = 0.0_dp 413 DO ispin = 1, nspin 414 DO ia = 1, natom 415 DO ib = ia + 1, natom 416 selnAB(ia, ib, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin) - occnumAB(ia, ib, ispin) 417 selnAB(ib, ia, ispin) = selnAB(ia, ib, ispin) 418 END DO 419 END DO 420 END DO 421 422 IF (.NOT. neglect_abc) THEN 423 ! calculate N_ABC 424 nabc = (natom*(natom - 1)*(natom - 2))/6 425 ALLOCATE (occnumABC(nabc, nspin)) 426 occnumABC = -1.0_dp 427 DO ispin = 1, nspin 428 CALL dbcsr_create(qmat_diag, name="MAO diagonal density", template=mao_dmat(1)%matrix) 429 CALL dbcsr_create(smat_diag, name="MAO diagonal overlap", template=mao_dmat(1)%matrix) 430 ! replicate the diagonal blocks of the density and overlap matrices 431 CALL dbcsr_get_block_diag(mao_qmat(ispin)%matrix, qmat_diag) 432 CALL dbcsr_replicate_all(qmat_diag) 433 CALL dbcsr_get_block_diag(mao_smat(ispin)%matrix, smat_diag) 434 CALL dbcsr_replicate_all(smat_diag) 435 iabc = 0 436 DO ia = 1, natom 437 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ia, col=ia, block=qblka, found=fo) 438 CPASSERT(fo) 439 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=sblka, found=fo) 440 CPASSERT(fo) 441 na = SIZE(qblka, 1) 442 DO ib = ia + 1, natom 443 ! screen with SEN(AB) 444 IF (selnAB(ia, ib, ispin) < eps_abc) THEN 445 iabc = iabc + (natom - ib) 446 CYCLE 447 END IF 448 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ib, col=ib, block=qblkb, found=fo) 449 CPASSERT(fo) 450 CALL dbcsr_get_block_p(matrix=smat_diag, row=ib, col=ib, block=sblkb, found=fo) 451 CPASSERT(fo) 452 nb = SIZE(qblkb, 1) 453 nab = na + nb 454 ALLOCATE (qmatab(na, nb), smatab(na, nb)) 455 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ib, & 456 block=block, found=found) 457 qmatab = 0.0_dp 458 IF (found) qmatab(1:na, 1:nb) = block(1:na, 1:nb) 459 CALL mp_sum(qmatab, para_env%group) 460 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ib, & 461 block=block, found=found) 462 smatab = 0.0_dp 463 IF (found) smatab(1:na, 1:nb) = block(1:na, 1:nb) 464 CALL mp_sum(smatab, para_env%group) 465 DO ic = ib + 1, natom 466 ! screen with SEN(AB) 467 IF ((selnAB(ia, ic, ispin) < eps_abc) .OR. (selnAB(ib, ic, ispin) < eps_abc)) THEN 468 iabc = iabc + 1 469 CYCLE 470 END IF 471 CALL dbcsr_get_block_p(matrix=qmat_diag, row=ic, col=ic, block=qblkc, found=fo) 472 CPASSERT(fo) 473 CALL dbcsr_get_block_p(matrix=smat_diag, row=ic, col=ic, block=sblkc, found=fo) 474 CPASSERT(fo) 475 nc = SIZE(qblkc, 1) 476 ALLOCATE (qmatac(na, nc), smatac(na, nc)) 477 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ia, col=ic, & 478 block=block, found=found) 479 qmatac = 0.0_dp 480 IF (found) qmatac(1:na, 1:nc) = block(1:na, 1:nc) 481 CALL mp_sum(qmatac, para_env%group) 482 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ia, col=ic, & 483 block=block, found=found) 484 smatac = 0.0_dp 485 IF (found) smatac(1:na, 1:nc) = block(1:na, 1:nc) 486 CALL mp_sum(smatac, para_env%group) 487 ALLOCATE (qmatbc(nb, nc), smatbc(nb, nc)) 488 CALL dbcsr_get_block_p(matrix=mao_qmat(ispin)%matrix, row=ib, col=ic, & 489 block=block, found=found) 490 qmatbc = 0.0_dp 491 IF (found) qmatbc(1:nb, 1:nc) = block(1:nb, 1:nc) 492 CALL mp_sum(qmatbc, para_env%group) 493 CALL dbcsr_get_block_p(matrix=mao_smat(ispin)%matrix, row=ib, col=ic, & 494 block=block, found=found) 495 smatbc = 0.0_dp 496 IF (found) smatbc(1:nb, 1:nc) = block(1:nb, 1:nc) 497 CALL mp_sum(smatbc, para_env%group) 498 ! 499 nabc = na + nb + nc 500 ALLOCATE (sab(nabc, nabc), sinv(nabc, nabc), qab(nabc, nabc)) 501 ! 502 qab(1:na, 1:na) = qblka(1:na, 1:na) 503 qab(na + 1:nab, na + 1:nab) = qblkb(1:nb, 1:nb) 504 qab(nab + 1:nabc, nab + 1:nabc) = qblkc(1:nc, 1:nc) 505 qab(1:na, na + 1:nab) = qmatab(1:na, 1:nb) 506 qab(na + 1:nab, 1:na) = TRANSPOSE(qmatab(1:na, 1:nb)) 507 qab(1:na, nab + 1:nabc) = qmatac(1:na, 1:nc) 508 qab(nab + 1:nabc, 1:na) = TRANSPOSE(qmatac(1:na, 1:nc)) 509 qab(na + 1:nab, nab + 1:nabc) = qmatbc(1:nb, 1:nc) 510 qab(nab + 1:nabc, na + 1:nab) = TRANSPOSE(qmatbc(1:nb, 1:nc)) 511 ! 512 sab(1:na, 1:na) = sblka(1:na, 1:na) 513 sab(na + 1:nab, na + 1:nab) = sblkb(1:nb, 1:nb) 514 sab(nab + 1:nabc, nab + 1:nabc) = sblkc(1:nc, 1:nc) 515 sab(1:na, na + 1:nab) = smatab(1:na, 1:nb) 516 sab(na + 1:nab, 1:na) = TRANSPOSE(smatab(1:na, 1:nb)) 517 sab(1:na, nab + 1:nabc) = smatac(1:na, 1:nc) 518 sab(nab + 1:nabc, 1:na) = TRANSPOSE(smatac(1:na, 1:nc)) 519 sab(na + 1:nab, nab + 1:nabc) = smatbc(1:nb, 1:nc) 520 sab(nab + 1:nabc, na + 1:nab) = TRANSPOSE(smatbc(1:nb, 1:nc)) 521 ! inv smat 522 sinv(1:nabc, 1:nabc) = sab(1:nabc, 1:nabc) 523 CALL invmat_symm(sinv) 524 ! Tr(Q*Sinv) 525 iabc = iabc + 1 526 me = MOD(iabc, para_env%num_pe) 527 IF (me == para_env%mepos) THEN 528 occnumABC(iabc, ispin) = SUM(qab*sinv) 529 ELSE 530 occnumABC(iabc, ispin) = 0.0_dp 531 END IF 532 ! 533 DEALLOCATE (sab, sinv, qab) 534 DEALLOCATE (qmatac, smatac) 535 DEALLOCATE (qmatbc, smatbc) 536 END DO 537 DEALLOCATE (qmatab, smatab) 538 END DO 539 END DO 540 CALL dbcsr_release(qmat_diag) 541 CALL dbcsr_release(smat_diag) 542 END DO 543 CALL mp_sum(occnumABC, para_env%group) 544 END IF 545 546 IF (.NOT. neglect_abc) THEN 547 ! calculate shared electron numbers (ABC) 548 nabc = (natom*(natom - 1)*(natom - 2))/6 549 ALLOCATE (selnABC(nabc, nspin)) 550 selnABC = 0.0_dp 551 DO ispin = 1, nspin 552 iabc = 0 553 DO ia = 1, natom 554 DO ib = ia + 1, natom 555 DO ic = ib + 1, natom 556 iabc = iabc + 1 557 IF (occnumABC(iabc, ispin) >= 0.0_dp) THEN 558 selnABC(iabc, ispin) = occnumA(ia, ispin) + occnumA(ib, ispin) + occnumA(ic, ispin) - & 559 occnumAB(ia, ib, ispin) - occnumAB(ia, ic, ispin) - occnumAB(ib, ic, ispin) + & 560 occnumABC(iabc, ispin) 561 END IF 562 END DO 563 END DO 564 END DO 565 END DO 566 END IF 567 568 ! calculate atomic charge 569 ALLOCATE (raq(natom, nspin)) 570 raq = 0.0_dp 571 DO ispin = 1, nspin 572 DO ia = 1, natom 573 raq(ia, ispin) = occnumA(ia, ispin) 574 DO ib = 1, natom 575 raq(ia, ispin) = raq(ia, ispin) - 0.5_dp*selnAB(ia, ib, ispin) 576 END DO 577 END DO 578 IF (.NOT. neglect_abc) THEN 579 iabc = 0 580 DO ia = 1, natom 581 DO ib = ia + 1, natom 582 DO ic = ib + 1, natom 583 iabc = iabc + 1 584 raq(ia, ispin) = raq(ia, ispin) + selnABC(iabc, ispin)/3._dp 585 raq(ib, ispin) = raq(ib, ispin) + selnABC(iabc, ispin)/3._dp 586 raq(ic, ispin) = raq(ic, ispin) + selnABC(iabc, ispin)/3._dp 587 END DO 588 END DO 589 END DO 590 END IF 591 END DO 592 593 ! calculate unassigned charge (from sum over atomic charges) 594 DO ispin = 1, nspin 595 deltaq = (electra(ispin) - SUM(raq(1:natom, ispin))) - ua_charge(ispin) 596 IF (unit_nr > 0) THEN 597 WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") & 598 "Cutoff error on charge", "Spin ", ispin, "error charge =", deltaq 599 END IF 600 END DO 601 602 ! analyze unassigned charge 603 ALLOCATE (uaq(natom, nspin)) 604 uaq = 0.0_dp 605 IF (analyze_ua) THEN 606 CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env) 607 CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, sab_all=sab_all) 608 CALL dbcsr_get_info(mao_coef(1)%matrix, row_blk_size=mao_blk_sizes, & 609 col_blk_size=col_blk_sizes, distribution=dbcsr_dist) 610 CALL dbcsr_get_info(matrix_s(1, 1)%matrix, row_blk_size=row_blk_sizes) 611 CALL dbcsr_create(amat, name="temp", template=matrix_s(1, 1)%matrix) 612 CALL dbcsr_create(tmat, name="temp", template=mao_coef(1)%matrix) 613 ! replicate diagonal of smm matrix 614 CALL dbcsr_get_block_diag(matrix_smm(1)%matrix, smat_diag) 615 CALL dbcsr_replicate_all(smat_diag) 616 617 ALLOCATE (orb_blk(natom), mao_blk(natom)) 618 DO ia = 1, natom 619 orb_blk = row_blk_sizes 620 mao_blk = row_blk_sizes 621 mao_blk(ia) = col_blk_sizes(ia) 622 CALL dbcsr_create(sumat, name="Smat", dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, & 623 row_blk_size=mao_blk, col_blk_size=mao_blk, nze=0) 624 CALL cp_dbcsr_alloc_block_from_nbl(sumat, sab_orb) 625 CALL dbcsr_create(cholmat, name="Cholesky matrix", dist=dbcsr_dist, & 626 matrix_type=dbcsr_type_no_symmetry, row_blk_size=mao_blk, col_blk_size=mao_blk, nze=0) 627 CALL dbcsr_create(rumat, name="Rmat", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, & 628 row_blk_size=orb_blk, col_blk_size=mao_blk, nze=0) 629 CALL cp_dbcsr_alloc_block_from_nbl(rumat, sab_orb, .TRUE.) 630 CALL dbcsr_create(crumat, name="Rmat*Umat", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, & 631 row_blk_size=orb_blk, col_blk_size=mao_blk, nze=0) 632 ! replicate row and col of smo matrix 633 ALLOCATE (rowblock(natom)) 634 DO ib = 1, natom 635 na = mao_blk_sizes(ia) 636 nb = row_blk_sizes(ib) 637 ALLOCATE (rowblock(ib)%mat(na, nb)) 638 rowblock(ib)%mat = 0.0_dp 639 CALL dbcsr_get_block_p(matrix=matrix_smo(1)%matrix, row=ia, col=ib, & 640 block=block, found=found) 641 IF (found) rowblock(ib)%mat(1:na, 1:nb) = block(1:na, 1:nb) 642 CALL mp_sum(rowblock(ib)%mat, para_env%group) 643 END DO 644 ! 645 DO ispin = 1, nspin 646 CALL dbcsr_copy(tmat, mao_coef(ispin)%matrix) 647 CALL dbcsr_replicate_all(tmat) 648 CALL dbcsr_iterator_start(dbcsr_iter, matrix_s(1, 1)%matrix) 649 DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter)) 650 CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, block) 651 CALL dbcsr_get_block_p(matrix=sumat, row=iatom, col=jatom, block=sblk, found=fos) 652 CPASSERT(fos) 653 CALL dbcsr_get_block_p(matrix=rumat, row=iatom, col=jatom, block=rblku, found=for) 654 CPASSERT(for) 655 CALL dbcsr_get_block_p(matrix=rumat, row=jatom, col=iatom, block=rblkl, found=for) 656 CPASSERT(for) 657 CALL dbcsr_get_block_p(matrix=tmat, row=ia, col=ia, block=cmao, found=found) 658 CPASSERT(found) 659 IF (iatom /= ia .AND. jatom /= ia) THEN 660 ! copy original overlap matrix 661 sblk = block 662 rblku = block 663 rblkl = TRANSPOSE(block) 664 ELSE IF (iatom /= ia) THEN 665 rblkl = TRANSPOSE(block) 666 sblk = MATMUL(TRANSPOSE(rowblock(iatom)%mat), cmao) 667 rblku = sblk 668 ELSE IF (jatom /= ia) THEN 669 rblku = block 670 sblk = MATMUL(TRANSPOSE(cmao), rowblock(jatom)%mat) 671 rblkl = TRANSPOSE(sblk) 672 ELSE 673 CALL dbcsr_get_block_p(matrix=smat_diag, row=ia, col=ia, block=block, found=found) 674 CPASSERT(found) 675 sblk = MATMUL(TRANSPOSE(cmao), MATMUL(block, cmao)) 676 rblku = MATMUL(TRANSPOSE(rowblock(ia)%mat), cmao) 677 END IF 678 END DO 679 CALL dbcsr_iterator_stop(dbcsr_iter) 680 ! Cholesky decomposition of SUMAT = U'U 681 CALL dbcsr_desymmetrize(sumat, cholmat) 682 CALL cp_dbcsr_cholesky_decompose(cholmat, para_env=para_env, blacs_env=blacs_env) 683 ! T = R*inv(U) 684 ssize = SUM(mao_blk) 685 CALL cp_dbcsr_cholesky_restore(rumat, ssize, cholmat, crumat, op="SOLVE", pos="RIGHT", & 686 transa="N", para_env=para_env, blacs_env=blacs_env) 687 ! A = T*transpose(T) 688 CALL dbcsr_multiply("N", "T", 1.0_dp, crumat, crumat, 0.0_dp, amat, & 689 filter_eps=eps_filter) 690 ! Tr(P*A) 691 CALL dbcsr_dot(matrix_p(ispin, 1)%matrix, amat, uaq(ia, ispin)) 692 uaq(ia, ispin) = uaq(ia, ispin) - electra(ispin) 693 END DO 694 ! 695 CALL dbcsr_release(sumat) 696 CALL dbcsr_release(cholmat) 697 CALL dbcsr_release(rumat) 698 CALL dbcsr_release(crumat) 699 ! 700 DO ib = 1, natom 701 DEALLOCATE (rowblock(ib)%mat) 702 END DO 703 DEALLOCATE (rowblock) 704 END DO 705 CALL dbcsr_release(smat_diag) 706 CALL dbcsr_release(amat) 707 CALL dbcsr_release(tmat) 708 DEALLOCATE (orb_blk, mao_blk) 709 END IF 710 ! 711 raq(1:natom, 1:nspin) = raq(1:natom, 1:nspin) - uaq(1:natom, 1:nspin) 712 DO ispin = 1, nspin 713 deltaq = electra(ispin) - SUM(raq(1:natom, ispin)) 714 IF (unit_nr > 0) THEN 715 WRITE (UNIT=unit_nr, FMT="(T2,A,T32,A,i2,T55,A,F12.8)") & 716 "Charge/Atom redistributed", "Spin ", ispin, "delta charge =", & 717 (deltaq + ua_charge(ispin))/REAL(natom, KIND=dp) 718 END IF 719 END DO 720 721 ! output charges 722 IF (unit_nr > 0) THEN 723 IF (nspin == 1) THEN 724 WRITE (unit_nr, "(/,T2,A,T40,A,T75,A)") "MAO atomic charges ", "Atom", "Charge" 725 ELSE 726 WRITE (unit_nr, "(/,T2,A,T40,A,T55,A,T70,A)") "MAO atomic charges ", "Atom", "Charge", "Spin Charge" 727 END IF 728 DO ispin = 1, nspin 729 deltaq = electra(ispin) - SUM(raq(1:natom, ispin)) 730 raq(:, ispin) = raq(:, ispin) + deltaq/REAL(natom, KIND=dp) 731 END DO 732 total_charge = 0.0_dp 733 total_spin = 0.0_dp 734 DO iatom = 1, natom 735 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 736 element_symbol=element_symbol, kind_number=ikind) 737 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff) 738 IF (nspin == 1) THEN 739 WRITE (unit_nr, "(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, zeff - raq(iatom, 1) 740 total_charge = total_charge + (zeff - raq(iatom, 1)) 741 ELSE 742 WRITE (unit_nr, "(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, & 743 zeff - raq(iatom, 1) - raq(iatom, 2), raq(iatom, 1) - raq(iatom, 2) 744 total_charge = total_charge + (zeff - raq(iatom, 1) - raq(iatom, 2)) 745 total_spin = total_spin + (raq(iatom, 1) - raq(iatom, 2)) 746 END IF 747 END DO 748 IF (nspin == 1) THEN 749 WRITE (unit_nr, "(T2,A,T69,F12.6)") "Total Charge", total_charge 750 ELSE 751 WRITE (unit_nr, "(T2,A,T49,F12.6,T69,F12.6)") "Total Charge", total_charge, total_spin 752 END IF 753 END IF 754 755 IF (analyze_ua) THEN 756 ! output unassigned charges 757 IF (unit_nr > 0) THEN 758 IF (nspin == 1) THEN 759 WRITE (unit_nr, "(/,T2,A,T40,A,T75,A)") "MAO hypervalent charges ", "Atom", "Charge" 760 ELSE 761 WRITE (unit_nr, "(/,T2,A,T40,A,T55,A,T70,A)") "MAO hypervalent charges ", "Atom", & 762 "Charge", "Spin Charge" 763 END IF 764 total_charge = 0.0_dp 765 total_spin = 0.0_dp 766 DO iatom = 1, natom 767 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 768 element_symbol=element_symbol) 769 IF (nspin == 1) THEN 770 WRITE (unit_nr, "(T30,I6,T42,A2,T69,F12.6)") iatom, element_symbol, uaq(iatom, 1) 771 total_charge = total_charge + uaq(iatom, 1) 772 ELSE 773 WRITE (unit_nr, "(T30,I6,T42,A2,T48,F12.6,T69,F12.6)") iatom, element_symbol, & 774 uaq(iatom, 1) + uaq(iatom, 2), uaq(iatom, 1) - uaq(iatom, 2) 775 total_charge = total_charge + uaq(iatom, 1) + uaq(iatom, 2) 776 total_spin = total_spin + uaq(iatom, 1) - uaq(iatom, 2) 777 END IF 778 END DO 779 IF (nspin == 1) THEN 780 WRITE (unit_nr, "(T2,A,T69,F12.6)") "Total Charge", total_charge 781 ELSE 782 WRITE (unit_nr, "(T2,A,T49,F12.6,T69,F12.6)") "Total Charge", total_charge, total_spin 783 END IF 784 END IF 785 END IF 786 787 ! output shared electron numbers AB 788 IF (unit_nr > 0) THEN 789 IF (nspin == 1) THEN 790 WRITE (unit_nr, "(/,T2,A,T31,A,T40,A,T78,A)") "Shared electron numbers ", "Atom", "Atom", "SEN" 791 ELSE 792 WRITE (unit_nr, "(/,T2,A,T31,A,T40,A,T51,A,T63,A,T71,A)") "Shared electron numbers ", "Atom", "Atom", & 793 "SEN(1)", "SEN(2)", "SEN(total)" 794 END IF 795 DO ia = 1, natom 796 DO ib = ia + 1, natom 797 CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa) 798 CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb) 799 IF (nspin == 1) THEN 800 IF (selnAB(ia, ib, 1) > eps_ab) THEN 801 WRITE (unit_nr, "(T26,I6,' ',A2,T35,I6,' ',A2,T69,F12.6)") ia, esa, ib, esb, selnAB(ia, ib, 1) 802 END IF 803 ELSE 804 IF ((selnAB(ia, ib, 1) + selnAB(ia, ib, 2)) > eps_ab) THEN 805 WRITE (unit_nr, "(T26,I6,' ',A2,T35,I6,' ',A2,T45,3F12.6)") ia, esa, ib, esb, & 806 selnAB(ia, ib, 1), selnAB(ia, ib, 2), (selnAB(ia, ib, 1) + selnAB(ia, ib, 2)) 807 END IF 808 END IF 809 END DO 810 END DO 811 END IF 812 813 IF (.NOT. neglect_abc) THEN 814 ! output shared electron numbers ABC 815 IF (unit_nr > 0) THEN 816 WRITE (unit_nr, "(/,T2,A,T40,A,T49,A,T58,A,T78,A)") "Shared electron numbers ABC", & 817 "Atom", "Atom", "Atom", "SEN" 818 senmax = 0.0_dp 819 iabc = 0 820 DO ia = 1, natom 821 DO ib = ia + 1, natom 822 DO ic = ib + 1, natom 823 iabc = iabc + 1 824 senabc = SUM(selnABC(iabc, :)) 825 senmax = MAX(senmax, senabc) 826 IF (senabc > eps_abc) THEN 827 CALL get_atomic_kind(atomic_kind=particle_set(ia)%atomic_kind, element_symbol=esa) 828 CALL get_atomic_kind(atomic_kind=particle_set(ib)%atomic_kind, element_symbol=esb) 829 CALL get_atomic_kind(atomic_kind=particle_set(ic)%atomic_kind, element_symbol=esc) 830 WRITE (unit_nr, "(T35,I6,' ',A2,T44,I6,' ',A2,T53,I6,' ',A2,T69,F12.6)") & 831 ia, esa, ib, esb, ic, esc, senabc 832 END IF 833 END DO 834 END DO 835 END DO 836 WRITE (unit_nr, "(T2,A,T69,F12.6)") "Maximum SEN value calculated", senmax 837 END IF 838 END IF 839 840 IF (unit_nr > 0) THEN 841 WRITE (unit_nr, '(/,T2,A)') & 842 '!---------------------------END OF MAO ANALYSIS-------------------------------!' 843 END IF 844 845 ! Deallocate temporary arrays 846 DEALLOCATE (occnumA, occnumAB, selnAB, raq, uaq) 847 IF (.NOT. neglect_abc) THEN 848 DEALLOCATE (occnumABC, selnABC) 849 END IF 850 851 ! Deallocate the neighbor list structure 852 CALL release_neighbor_list_sets(smm_list) 853 CALL release_neighbor_list_sets(smo_list) 854 855 DEALLOCATE (mao_basis_set_list, orb_basis_set_list) 856 857 IF (ASSOCIATED(matrix_smm)) CALL dbcsr_deallocate_matrix_set(matrix_smm) 858 IF (ASSOCIATED(matrix_smo)) CALL dbcsr_deallocate_matrix_set(matrix_smo) 859 IF (ASSOCIATED(matrix_q)) CALL dbcsr_deallocate_matrix_set(matrix_q) 860 861 IF (ASSOCIATED(mao_coef)) CALL dbcsr_deallocate_matrix_set(mao_coef) 862 IF (ASSOCIATED(mao_dmat)) CALL dbcsr_deallocate_matrix_set(mao_dmat) 863 IF (ASSOCIATED(mao_smat)) CALL dbcsr_deallocate_matrix_set(mao_smat) 864 IF (ASSOCIATED(mao_qmat)) CALL dbcsr_deallocate_matrix_set(mao_qmat) 865 866 CALL timestop(handle) 867 868 END SUBROUTINE mao_analysis 869 870END MODULE mao_wfn_analysis 871