LCOV - code coverage report
Current view: top level - src - mao_wfn_analysis.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 98.4 % 510 502
Test Date: 2025-07-25 12:55:17 Functions: 33.3 % 3 1

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

Generated by: LCOV version 2.0-1