LCOV - code coverage report
Current view: top level - src - population_analyses.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 265 271 97.8 %
Date: 2024-03-28 07:31:50 Functions: 4 4 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : ! **************************************************************************************************
       8             : !> \brief Provide various population analyses and print the requested output
       9             : !>        information
      10             : !>
      11             : !> \author  Matthias Krack (MK)
      12             : !> \date    09.07.2010
      13             : !> \version 1.0
      14             : ! **************************************************************************************************
      15             : 
      16             : MODULE population_analyses
      17             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      18             :                                               get_atomic_kind,&
      19             :                                               get_atomic_kind_set
      20             :    USE basis_set_types,                 ONLY: get_gto_basis_set,&
      21             :                                               gto_basis_set_type
      22             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      23             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      24             :                                               cp_dbcsr_sm_fm_multiply
      25             :    USE cp_dbcsr_output,                 ONLY: cp_dbcsr_write_sparse_matrix,&
      26             :                                               write_fm_with_basis_info
      27             :    USE cp_fm_diag,                      ONLY: cp_fm_power
      28             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      29             :                                               cp_fm_struct_release,&
      30             :                                               cp_fm_struct_type
      31             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      32             :                                               cp_fm_get_diag,&
      33             :                                               cp_fm_release,&
      34             :                                               cp_fm_type
      35             :    USE cp_result_methods,               ONLY: cp_results_erase,&
      36             :                                               put_results
      37             :    USE cp_result_types,                 ONLY: cp_result_type
      38             :    USE dbcsr_api,                       ONLY: &
      39             :         dbcsr_copy, dbcsr_deallocate_matrix, dbcsr_get_block_p, dbcsr_iterator_blocks_left, &
      40             :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      41             :         dbcsr_p_type, dbcsr_set, dbcsr_setname, dbcsr_type
      42             :    USE kinds,                           ONLY: default_string_length,&
      43             :                                               dp
      44             :    USE machine,                         ONLY: m_flush
      45             :    USE message_passing,                 ONLY: mp_para_env_type
      46             :    USE orbital_pointers,                ONLY: nso
      47             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      48             :    USE particle_methods,                ONLY: get_particle_set
      49             :    USE particle_types,                  ONLY: particle_type
      50             :    USE qs_environment_types,            ONLY: get_qs_env,&
      51             :                                               qs_environment_type
      52             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      53             :                                               get_qs_kind_set,&
      54             :                                               qs_kind_type
      55             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      56             :                                               qs_rho_type
      57             :    USE scf_control_types,               ONLY: scf_control_type
      58             : #include "./base/base_uses.f90"
      59             : 
      60             :    IMPLICIT NONE
      61             : 
      62             :    PRIVATE
      63             : 
      64             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'population_analyses'
      65             : 
      66             :    PUBLIC :: lowdin_population_analysis, &
      67             :              mulliken_population_analysis
      68             : 
      69             : CONTAINS
      70             : 
      71             : ! **************************************************************************************************
      72             : !> \brief Perform a Lowdin population analysis based on a symmetric
      73             : !>        orthogonalisation of the density matrix using S^(1/2)
      74             : !>
      75             : !> \param qs_env ...
      76             : !> \param output_unit ...
      77             : !> \param print_level ...
      78             : !> \date    06.07.2010
      79             : !> \author  Matthias Krack (MK)
      80             : !> \version 1.0
      81             : ! **************************************************************************************************
      82          82 :    SUBROUTINE lowdin_population_analysis(qs_env, output_unit, print_level)
      83             : 
      84             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      85             :       INTEGER, INTENT(IN)                                :: output_unit, print_level
      86             : 
      87             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'lowdin_population_analysis'
      88             : 
      89             :       CHARACTER(LEN=default_string_length)               :: headline
      90             :       INTEGER                                            :: handle, ispin, ndep, nsgf, nspin
      91             :       LOGICAL                                            :: print_gop
      92          82 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: orbpop
      93          82 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      94             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
      95             :       TYPE(cp_fm_struct_type), POINTER                   :: fmstruct
      96             :       TYPE(cp_fm_type)                                   :: fm_s_half, fm_work1, fm_work2
      97          82 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p
      98          82 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrixkp_p, matrixkp_s
      99             :       TYPE(dbcsr_type), POINTER                          :: sm_p, sm_s
     100             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     101          82 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     102          82 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     103             :       TYPE(qs_rho_type), POINTER                         :: rho
     104             :       TYPE(scf_control_type), POINTER                    :: scf_control
     105             : 
     106          82 :       CALL timeset(routineN, handle)
     107             : 
     108          82 :       NULLIFY (atomic_kind_set)
     109          82 :       NULLIFY (qs_kind_set)
     110          82 :       NULLIFY (fmstruct)
     111          82 :       NULLIFY (matrix_p)
     112          82 :       NULLIFY (matrixkp_p)
     113          82 :       NULLIFY (matrixkp_s)
     114          82 :       NULLIFY (orbpop)
     115          82 :       NULLIFY (particle_set)
     116          82 :       NULLIFY (rho)
     117          82 :       NULLIFY (scf_control)
     118          82 :       NULLIFY (sm_p)
     119          82 :       NULLIFY (sm_s)
     120             :       NULLIFY (orbpop)
     121          82 :       NULLIFY (para_env)
     122          82 :       NULLIFY (blacs_env)
     123             : 
     124             :       CALL get_qs_env(qs_env=qs_env, &
     125             :                       atomic_kind_set=atomic_kind_set, &
     126             :                       qs_kind_set=qs_kind_set, &
     127             :                       matrix_s_kp=matrixkp_s, &
     128             :                       particle_set=particle_set, &
     129             :                       rho=rho, &
     130             :                       scf_control=scf_control, &
     131             :                       para_env=para_env, &
     132          82 :                       blacs_env=blacs_env)
     133             : 
     134          82 :       CPASSERT(ASSOCIATED(atomic_kind_set))
     135          82 :       CPASSERT(ASSOCIATED(qs_kind_set))
     136          82 :       CPASSERT(ASSOCIATED(matrixkp_s))
     137          82 :       CPASSERT(ASSOCIATED(particle_set))
     138          82 :       CPASSERT(ASSOCIATED(rho))
     139          82 :       CPASSERT(ASSOCIATED(scf_control))
     140             : 
     141          82 :       IF (SIZE(matrixkp_s, 2) > 1) THEN
     142             : 
     143           0 :          CPWARN("Lowdin population analysis not implemented for k-points.")
     144             : 
     145             :       ELSE
     146             : 
     147          82 :          sm_s => matrixkp_s(1, 1)%matrix ! Overlap matrix in sparse format
     148          82 :          CALL qs_rho_get(rho, rho_ao_kp=matrixkp_p) ! Density matrices in sparse format
     149             : 
     150          82 :          matrix_p => matrixkp_p(:, 1)
     151          82 :          nspin = SIZE(matrix_p, 1)
     152             : 
     153             :          ! Get the total number of contracted spherical Gaussian basis functions
     154          82 :          CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
     155             : 
     156             :          ! Provide an array to store the orbital populations for each spin
     157         328 :          ALLOCATE (orbpop(nsgf, nspin))
     158        2612 :          orbpop(:, :) = 0.0_dp
     159             : 
     160             :          ! Write headline
     161          82 :          IF (output_unit > 0) THEN
     162          41 :             WRITE (UNIT=output_unit, FMT="(/,/,T2,A)") "LOWDIN POPULATION ANALYSIS"
     163             :          END IF
     164             : 
     165             :          ! Provide full size work matrices
     166             :          CALL cp_fm_struct_create(fmstruct=fmstruct, &
     167             :                                   para_env=para_env, &
     168             :                                   context=blacs_env, &
     169             :                                   nrow_global=nsgf, &
     170          82 :                                   ncol_global=nsgf)
     171             :          CALL cp_fm_create(matrix=fm_s_half, &
     172             :                            matrix_struct=fmstruct, &
     173          82 :                            name="S^(1/2) MATRIX")
     174             :          CALL cp_fm_create(matrix=fm_work1, &
     175             :                            matrix_struct=fmstruct, &
     176          82 :                            name="FULL WORK MATRIX 1")
     177          82 :          headline = "SYMMETRICALLY ORTHOGONALISED DENSITY MATRIX"
     178             :          CALL cp_fm_create(matrix=fm_work2, &
     179             :                            matrix_struct=fmstruct, &
     180          82 :                            name=TRIM(headline))
     181          82 :          CALL cp_fm_struct_release(fmstruct=fmstruct)
     182             : 
     183             :          ! Build full S^(1/2) matrix (computationally expensive)
     184          82 :          CALL copy_dbcsr_to_fm(sm_s, fm_s_half)
     185          82 :          CALL cp_fm_power(fm_s_half, fm_work1, 0.5_dp, scf_control%eps_eigval, ndep)
     186          82 :          IF (ndep /= 0) &
     187             :             CALL cp_warn(__LOCATION__, &
     188             :                          "Overlap matrix exhibits linear dependencies. At least some "// &
     189           0 :                          "eigenvalues have been quenched.")
     190             : 
     191             :          ! Build Lowdin population matrix for each spin
     192         168 :          DO ispin = 1, nspin
     193          86 :             sm_p => matrix_p(ispin)%matrix ! Density matrix for spin ispin in sparse format
     194             :             ! Calculate S^(1/2)*P*S^(1/2) as a full matrix (Lowdin)
     195          86 :             CALL cp_dbcsr_sm_fm_multiply(sm_p, fm_s_half, fm_work1, nsgf)
     196             :             CALL parallel_gemm(transa="N", &
     197             :                                transb="N", &
     198             :                                m=nsgf, &
     199             :                                n=nsgf, &
     200             :                                k=nsgf, &
     201             :                                alpha=1.0_dp, &
     202             :                                matrix_a=fm_s_half, &
     203             :                                matrix_b=fm_work1, &
     204             :                                beta=0.0_dp, &
     205          86 :                                matrix_c=fm_work2)
     206          86 :             IF (print_level > 2) THEN
     207             :                ! Write the full Lowdin population matrix
     208           4 :                IF (nspin > 1) THEN
     209           4 :                   IF (ispin == 1) THEN
     210           2 :                      fm_work2%name = TRIM(headline)//" FOR ALPHA SPIN"
     211             :                   ELSE
     212           2 :                      fm_work2%name = TRIM(headline)//" FOR BETA SPIN"
     213             :                   END IF
     214             :                END IF
     215             :                CALL write_fm_with_basis_info(fm_work2, 4, 6, qs_env, para_env, &
     216           4 :                                              output_unit=output_unit)
     217             :             END IF
     218         168 :             CALL cp_fm_get_diag(fm_work2, orbpop(:, ispin))
     219             :          END DO ! next spin ispin
     220             : 
     221             :          ! Write atomic populations and charges
     222          82 :          IF (output_unit > 0) THEN
     223          41 :             print_gop = (print_level > 1) ! Print also orbital populations
     224          41 :             CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
     225             :          END IF
     226             : 
     227             :          ! Release local working storage
     228          82 :          CALL cp_fm_release(matrix=fm_s_half)
     229          82 :          CALL cp_fm_release(matrix=fm_work1)
     230          82 :          CALL cp_fm_release(matrix=fm_work2)
     231         164 :          IF (ASSOCIATED(orbpop)) THEN
     232          82 :             DEALLOCATE (orbpop)
     233             :          END IF
     234             : 
     235             :       END IF
     236             : 
     237          82 :       CALL timestop(handle)
     238             : 
     239          82 :    END SUBROUTINE lowdin_population_analysis
     240             : 
     241             : ! **************************************************************************************************
     242             : !> \brief Perform a Mulliken population analysis
     243             : !>
     244             : !> \param qs_env ...
     245             : !> \param output_unit ...
     246             : !> \param print_level ...
     247             : !> \date    10.07.2010
     248             : !> \author  Matthias Krack (MK)
     249             : !> \version 1.0
     250             : ! **************************************************************************************************
     251        4488 :    SUBROUTINE mulliken_population_analysis(qs_env, output_unit, print_level)
     252             : 
     253             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     254             :       INTEGER, INTENT(IN)                                :: output_unit, print_level
     255             : 
     256             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'mulliken_population_analysis'
     257             : 
     258             :       CHARACTER(LEN=default_string_length)               :: headline
     259             :       INTEGER                                            :: blk, handle, iatom, ic, isgf, ispin, &
     260             :                                                             jatom, jsgf, natom, nsgf, nspin, sgfa, &
     261             :                                                             sgfb
     262             :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: first_sgf_atom
     263             :       LOGICAL                                            :: found, print_gop
     264             :       REAL(KIND=dp)                                      :: ps
     265        4488 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: orbpop, p_block, ps_block, s_block
     266        4488 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     267             :       TYPE(dbcsr_iterator_type)                          :: iter
     268        4488 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_p, matrix_s
     269             :       TYPE(dbcsr_type), POINTER                          :: sm_p, sm_ps, sm_s
     270             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     271        4488 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     272        4488 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     273             :       TYPE(qs_rho_type), POINTER                         :: rho
     274             : 
     275        4488 :       CALL timeset(routineN, handle)
     276             : 
     277        4488 :       NULLIFY (atomic_kind_set)
     278        4488 :       NULLIFY (qs_kind_set)
     279        4488 :       NULLIFY (matrix_p)
     280        4488 :       NULLIFY (matrix_s)
     281             :       NULLIFY (orbpop)
     282        4488 :       NULLIFY (particle_set)
     283        4488 :       NULLIFY (ps_block)
     284        4488 :       NULLIFY (p_block)
     285        4488 :       NULLIFY (rho)
     286        4488 :       NULLIFY (sm_p)
     287        4488 :       NULLIFY (sm_ps)
     288        4488 :       NULLIFY (sm_s)
     289        4488 :       NULLIFY (s_block)
     290        4488 :       NULLIFY (para_env)
     291             : 
     292             :       CALL get_qs_env(qs_env=qs_env, &
     293             :                       atomic_kind_set=atomic_kind_set, &
     294             :                       qs_kind_set=qs_kind_set, &
     295             :                       matrix_s_kp=matrix_s, &
     296             :                       particle_set=particle_set, &
     297             :                       rho=rho, &
     298        4488 :                       para_env=para_env)
     299             : 
     300        4488 :       CPASSERT(ASSOCIATED(atomic_kind_set))
     301        4488 :       CPASSERT(ASSOCIATED(qs_kind_set))
     302        4488 :       CPASSERT(ASSOCIATED(particle_set))
     303        4488 :       CPASSERT(ASSOCIATED(rho))
     304        4488 :       CPASSERT(ASSOCIATED(matrix_s))
     305             : 
     306        4488 :       CALL qs_rho_get(rho, rho_ao_kp=matrix_p) ! Density matrices in sparse format
     307        4488 :       nspin = SIZE(matrix_p, 1)
     308             : 
     309             :       ! Get the total number of contracted spherical Gaussian basis functions
     310        4488 :       CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
     311        4488 :       CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
     312       13464 :       ALLOCATE (first_sgf_atom(natom))
     313       23848 :       first_sgf_atom(:) = 0
     314             : 
     315        4488 :       CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf_atom)
     316             : 
     317             :       ! Provide an array to store the orbital populations for each spin
     318       17952 :       ALLOCATE (orbpop(nsgf, nspin))
     319      162030 :       orbpop(:, :) = 0.0_dp
     320             : 
     321             :       ! Write headline
     322        4488 :       IF (output_unit > 0) THEN
     323             :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     324        2258 :             '!-----------------------------------------------------------------------------!'
     325        2258 :          WRITE (UNIT=output_unit, FMT="(T22,A)") "Mulliken Population Analysis"
     326             :       END IF
     327             : 
     328             :       ! Create a DBCSR work matrix, if needed
     329        4488 :       IF (print_level > 2) THEN
     330           2 :          sm_s => matrix_s(1, 1)%matrix ! Overlap matrix in sparse format
     331           2 :          ALLOCATE (sm_ps)
     332           2 :          headline = "MULLIKEN NET ATOMIC ORBITAL AND OVERLAP POPULATION MATRIX"
     333             :          CALL dbcsr_copy(matrix_b=sm_ps, &
     334             :                          matrix_a=sm_s, &
     335           2 :                          name=TRIM(headline))
     336             :       END IF
     337             : 
     338             :       ! Build Mulliken population matrix for each spin
     339        9720 :       DO ispin = 1, nspin
     340       14518 :          DO ic = 1, SIZE(matrix_s, 2)
     341        9286 :             IF (print_level > 2) THEN
     342           4 :                CALL dbcsr_set(sm_ps, 0.0_dp)
     343             :             END IF
     344        9286 :             sm_s => matrix_s(1, ic)%matrix ! Overlap matrix in sparse format
     345        9286 :             sm_p => matrix_p(ispin, ic)%matrix ! Density matrix for spin ispin in sparse format
     346             :             ! Calculate Hadamard product of P and S as sparse matrix (Mulliken)
     347             :             ! CALL dbcsr_hadamard_product(sm_p,sm_s,sm_ps)
     348        9286 :             CALL dbcsr_iterator_start(iter, sm_s)
     349       94597 :             DO WHILE (dbcsr_iterator_blocks_left(iter))
     350       85311 :                CALL dbcsr_iterator_next_block(iter, iatom, jatom, s_block, blk)
     351       85311 :                IF (.NOT. (ASSOCIATED(s_block))) CYCLE
     352             :                CALL dbcsr_get_block_p(matrix=sm_p, &
     353             :                                       row=iatom, &
     354             :                                       col=jatom, &
     355             :                                       block=p_block, &
     356       85311 :                                       found=found)
     357       85311 :                IF (print_level > 2) THEN
     358             :                   CALL dbcsr_get_block_p(matrix=sm_ps, &
     359             :                                          row=iatom, &
     360             :                                          col=jatom, &
     361             :                                          block=ps_block, &
     362          12 :                                          found=found)
     363          12 :                   CPASSERT(ASSOCIATED(ps_block))
     364             :                END IF
     365             : 
     366       85311 :                sgfb = first_sgf_atom(jatom)
     367      462017 :                DO jsgf = 1, SIZE(s_block, 2)
     368     4842278 :                   DO isgf = 1, SIZE(s_block, 1)
     369     4465572 :                      ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
     370     4465572 :                      IF (ASSOCIATED(ps_block)) ps_block(isgf, jsgf) = ps_block(isgf, jsgf) + ps
     371     4842278 :                      orbpop(sgfb, ispin) = orbpop(sgfb, ispin) + ps
     372             :                   END DO
     373      462017 :                   sgfb = sgfb + 1
     374             :                END DO
     375       94597 :                IF (iatom /= jatom) THEN
     376       69504 :                   sgfa = first_sgf_atom(iatom)
     377      364161 :                   DO isgf = 1, SIZE(s_block, 1)
     378     3221243 :                      DO jsgf = 1, SIZE(s_block, 2)
     379     2926586 :                         ps = p_block(isgf, jsgf)*s_block(isgf, jsgf)
     380     3221243 :                         orbpop(sgfa, ispin) = orbpop(sgfa, ispin) + ps
     381             :                      END DO
     382      364161 :                      sgfa = sgfa + 1
     383             :                   END DO
     384             :                END IF
     385             :             END DO
     386       14518 :             CALL dbcsr_iterator_stop(iter)
     387             :          END DO
     388             : 
     389        9720 :          IF (print_level > 2) THEN
     390             :             ! Write the full Mulliken net AO and overlap population matrix
     391           4 :             IF (nspin > 1) THEN
     392           4 :                IF (ispin == 1) THEN
     393           2 :                   CALL dbcsr_setname(sm_ps, TRIM(headline)//" For Alpha Spin")
     394             :                ELSE
     395           2 :                   CALL dbcsr_setname(sm_ps, TRIM(headline)//" For Beta Spin")
     396             :                END IF
     397             :             END IF
     398             :             CALL cp_dbcsr_write_sparse_matrix(sm_ps, 4, 6, qs_env, para_env, &
     399           4 :                                               output_unit=output_unit)
     400             :          END IF
     401             :       END DO
     402             : 
     403      319572 :       CALL para_env%sum(orbpop)
     404             : 
     405             :       ! Write atomic populations and charges
     406        4488 :       IF (output_unit > 0) THEN
     407        2258 :          print_gop = (print_level > 1) ! Print also orbital populations
     408        2258 :          CALL write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, print_gop)
     409             :       END IF
     410             : 
     411             :       ! Save the Mulliken charges to results
     412        4488 :       CALL save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
     413             : 
     414             :       ! Release local working storage
     415        4488 :       IF (ASSOCIATED(sm_ps)) CALL dbcsr_deallocate_matrix(sm_ps)
     416        4488 :       IF (ASSOCIATED(orbpop)) THEN
     417        4488 :          DEALLOCATE (orbpop)
     418             :       END IF
     419        4488 :       IF (ALLOCATED(first_sgf_atom)) THEN
     420        4488 :          DEALLOCATE (first_sgf_atom)
     421             :       END IF
     422             : 
     423        4488 :       IF (output_unit > 0) THEN
     424             :          WRITE (UNIT=output_unit, FMT="(T2,A)") &
     425        2258 :             '!-----------------------------------------------------------------------------!'
     426             :       END IF
     427             : 
     428        4488 :       CALL timestop(handle)
     429             : 
     430       13464 :    END SUBROUTINE mulliken_population_analysis
     431             : 
     432             : ! **************************************************************************************************
     433             : !> \brief Save the Mulliken charges in results
     434             : !>
     435             : !> \param orbpop ...
     436             : !> \param atomic_kind_set ...
     437             : !> \param qs_kind_set ...
     438             : !> \param particle_set ...
     439             : !> \param qs_env ...
     440             : !> \date    27.05.2022
     441             : !> \author  Bo Thomsen (BT)
     442             : !> \version 1.0
     443             : ! **************************************************************************************************
     444        4488 :    SUBROUTINE save_mulliken_charges(orbpop, atomic_kind_set, qs_kind_set, particle_set, qs_env)
     445             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: orbpop
     446             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     447             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     448             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     449             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     450             : 
     451             :       CHARACTER(LEN=default_string_length)               :: description
     452             :       INTEGER                                            :: iao, iatom, ikind, iset, isgf, ishell, &
     453             :                                                             iso, l, natom, nset, nsgf, nspin
     454        4488 :       INTEGER, DIMENSION(:), POINTER                     :: nshell
     455        4488 :       INTEGER, DIMENSION(:, :), POINTER                  :: lshell
     456             :       REAL(KIND=dp)                                      :: zeff
     457             :       REAL(KIND=dp), DIMENSION(3)                        :: sumorbpop
     458             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: charges_save
     459             :       TYPE(cp_result_type), POINTER                      :: results
     460             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     461             : 
     462        4488 :       NULLIFY (lshell)
     463        4488 :       NULLIFY (nshell)
     464        4488 :       NULLIFY (orb_basis_set)
     465             : 
     466           0 :       CPASSERT(ASSOCIATED(orbpop))
     467        4488 :       CPASSERT(ASSOCIATED(atomic_kind_set))
     468        4488 :       CPASSERT(ASSOCIATED(particle_set))
     469             : 
     470        4488 :       nspin = SIZE(orbpop, 2)
     471             : 
     472        4488 :       CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
     473        4488 :       CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
     474        4488 :       CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
     475        4488 :       NULLIFY (results)
     476        4488 :       CALL get_qs_env(qs_env, results=results)
     477       13464 :       ALLOCATE (charges_save(natom))
     478             : 
     479        4488 :       iao = 1
     480       23848 :       DO iatom = 1, natom
     481       19360 :          sumorbpop(:) = 0.0_dp
     482       19360 :          NULLIFY (orb_basis_set)
     483             :          CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
     484       19360 :                               kind_number=ikind)
     485       19360 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
     486       23848 :          IF (ASSOCIATED(orb_basis_set)) THEN
     487             :             CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     488             :                                    nset=nset, &
     489             :                                    nshell=nshell, &
     490       19360 :                                    l=lshell)
     491       19360 :             isgf = 1
     492       52520 :             DO iset = 1, nset
     493      112722 :                DO ishell = 1, nshell(iset)
     494       60202 :                   l = lshell(ishell, iset)
     495      221928 :                   DO iso = 1, nso(l)
     496      128566 :                      IF (nspin == 1) THEN
     497      104822 :                         sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
     498             :                      ELSE
     499       71232 :                         sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
     500       23744 :                         sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
     501             :                      END IF
     502      128566 :                      isgf = isgf + 1
     503      188768 :                      iao = iao + 1
     504             :                   END DO
     505             :                END DO
     506             :             END DO
     507       19360 :             IF (nspin == 1) THEN
     508       16502 :                charges_save(iatom) = zeff - sumorbpop(1)
     509             :             ELSE
     510        2858 :                charges_save(iatom) = zeff - sumorbpop(1) - sumorbpop(2)
     511             :             END IF
     512             :          END IF ! atom has an orbital basis
     513             :       END DO ! next atom iatom
     514             : 
     515             :       ! Store charges in results
     516        4488 :       description = "[MULLIKEN-CHARGES]"
     517        4488 :       CALL cp_results_erase(results=results, description=description)
     518             :       CALL put_results(results=results, description=description, &
     519        4488 :                        values=charges_save)
     520             : 
     521        4488 :       DEALLOCATE (charges_save)
     522             : 
     523        8976 :    END SUBROUTINE save_mulliken_charges
     524             : 
     525             : ! **************************************************************************************************
     526             : !> \brief Write atomic orbital populations and net atomic charges
     527             : !>
     528             : !> \param orbpop ...
     529             : !> \param atomic_kind_set ...
     530             : !> \param qs_kind_set ...
     531             : !> \param particle_set ...
     532             : !> \param output_unit ...
     533             : !> \param print_orbital_contributions ...
     534             : !> \date    07.07.2010
     535             : !> \author  Matthias Krack (MK)
     536             : !> \version 1.0
     537             : ! **************************************************************************************************
     538        6897 :    SUBROUTINE write_orbpop(orbpop, atomic_kind_set, qs_kind_set, particle_set, output_unit, &
     539             :                            print_orbital_contributions)
     540             : 
     541             :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: orbpop
     542             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     543             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     544             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     545             :       INTEGER, INTENT(IN)                                :: output_unit
     546             :       LOGICAL, INTENT(IN)                                :: print_orbital_contributions
     547             : 
     548             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'write_orbpop'
     549             : 
     550             :       CHARACTER(LEN=2)                                   :: element_symbol
     551        2299 :       CHARACTER(LEN=6), DIMENSION(:), POINTER            :: sgf_symbol
     552             :       INTEGER                                            :: handle, iao, iatom, ikind, iset, isgf, &
     553             :                                                             ishell, iso, l, natom, nset, nsgf, &
     554             :                                                             nspin
     555        2299 :       INTEGER, DIMENSION(:), POINTER                     :: nshell
     556        2299 :       INTEGER, DIMENSION(:, :), POINTER                  :: lshell
     557             :       REAL(KIND=dp)                                      :: zeff
     558             :       REAL(KIND=dp), DIMENSION(3)                        :: sumorbpop, totsumorbpop
     559             :       TYPE(gto_basis_set_type), POINTER                  :: orb_basis_set
     560             : 
     561        2299 :       CALL timeset(routineN, handle)
     562             : 
     563        2299 :       NULLIFY (lshell)
     564        2299 :       NULLIFY (nshell)
     565        2299 :       NULLIFY (orb_basis_set)
     566        2299 :       NULLIFY (sgf_symbol)
     567             : 
     568        2299 :       CPASSERT(ASSOCIATED(orbpop))
     569        2299 :       CPASSERT(ASSOCIATED(atomic_kind_set))
     570        2299 :       CPASSERT(ASSOCIATED(particle_set))
     571             : 
     572        2299 :       nspin = SIZE(orbpop, 2)
     573             : 
     574        2299 :       CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
     575        2299 :       CALL get_qs_kind_set(qs_kind_set, nsgf=nsgf)
     576             : 
     577             :       ! Select and write headline
     578        2299 :       IF (nspin == 1) THEN
     579        1925 :          IF (print_orbital_contributions) THEN
     580             :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     581           0 :                "# Orbital  AO symbol  Orbital population                            Net charge"
     582             :          ELSE
     583             :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     584        1925 :                "#  Atom  Element  Kind  Atomic population                           Net charge"
     585             :          END IF
     586             :       ELSE
     587         374 :          IF (print_orbital_contributions) THEN
     588             :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     589           2 :                "# Orbital  AO symbol  Orbital population (alpha,beta)  Net charge  Spin moment"
     590             :          ELSE
     591             :             WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     592         372 :                "#  Atom  Element  Kind  Atomic population (alpha,beta) Net charge  Spin moment"
     593             :          END IF
     594             :       END IF
     595             : 
     596        2299 :       totsumorbpop(:) = 0.0_dp
     597             : 
     598        2299 :       iao = 1
     599       12118 :       DO iatom = 1, natom
     600        9819 :          sumorbpop(:) = 0.0_dp
     601        9819 :          NULLIFY (orb_basis_set)
     602             :          CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
     603             :                               element_symbol=element_symbol, &
     604        9819 :                               kind_number=ikind)
     605        9819 :          CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, zeff=zeff)
     606       12118 :          IF (ASSOCIATED(orb_basis_set)) THEN
     607             :             CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
     608             :                                    nset=nset, &
     609             :                                    nshell=nshell, &
     610             :                                    l=lshell, &
     611        9819 :                                    sgf_symbol=sgf_symbol)
     612        9819 :             isgf = 1
     613       26850 :             DO iset = 1, nset
     614       57532 :                DO ishell = 1, nshell(iset)
     615       30682 :                   l = lshell(ishell, iset)
     616      113293 :                   DO iso = 1, nso(l)
     617       65580 :                      IF (nspin == 1) THEN
     618       53643 :                         sumorbpop(1) = sumorbpop(1) + orbpop(iao, 1)
     619       53643 :                         IF (print_orbital_contributions) THEN
     620           0 :                            IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
     621             :                            WRITE (UNIT=output_unit, &
     622             :                                   FMT="(T2,I9,2X,A2,1X,A,T30,F12.6)") &
     623           0 :                               iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1)
     624             :                         END IF
     625             :                      ELSE
     626       35811 :                         sumorbpop(1:2) = sumorbpop(1:2) + orbpop(iao, 1:2)
     627       11937 :                         sumorbpop(3) = sumorbpop(3) + orbpop(iao, 1) - orbpop(iao, 2)
     628       11937 :                         IF (print_orbital_contributions) THEN
     629          78 :                            IF (isgf == 1) WRITE (UNIT=output_unit, FMT="(A)") ""
     630             :                            WRITE (UNIT=output_unit, &
     631             :                                   FMT="(T2,I9,2X,A2,1X,A,T29,2(1X,F12.6),T68,F12.6)") &
     632          78 :                               iao, element_symbol, sgf_symbol(isgf), orbpop(iao, 1:2), &
     633         156 :                               orbpop(iao, 1) - orbpop(iao, 2)
     634             :                         END IF
     635             :                      END IF
     636       65580 :                      isgf = isgf + 1
     637       96262 :                      iao = iao + 1
     638             :                   END DO
     639             :                END DO
     640             :             END DO
     641        9819 :             IF (nspin == 1) THEN
     642        8385 :                totsumorbpop(1) = totsumorbpop(1) + sumorbpop(1)
     643        8385 :                totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1)
     644             :                WRITE (UNIT=output_unit, &
     645             :                       FMT="(T2,I7,5X,A2,2X,I6,T30,F12.6,T68,F12.6)") &
     646        8385 :                   iatom, element_symbol, ikind, sumorbpop(1), zeff - sumorbpop(1)
     647             :             ELSE
     648        4302 :                totsumorbpop(1:2) = totsumorbpop(1:2) + sumorbpop(1:2)
     649        1434 :                totsumorbpop(3) = totsumorbpop(3) + zeff - sumorbpop(1) - sumorbpop(2)
     650             :                WRITE (UNIT=output_unit, &
     651             :                       FMT="(T2,I7,5X,A2,2X,I6,T28,4(1X,F12.6))") &
     652        1434 :                   iatom, element_symbol, ikind, sumorbpop(1:2), &
     653        2868 :                   zeff - sumorbpop(1) - sumorbpop(2), sumorbpop(3)
     654             :             END IF
     655             :          END IF ! atom has an orbital basis
     656             :       END DO ! next atom iatom
     657             : 
     658             :       ! Write total sums
     659        2299 :       IF (print_orbital_contributions) WRITE (UNIT=output_unit, FMT="(A)") ""
     660        2299 :       IF (nspin == 1) THEN
     661             :          WRITE (UNIT=output_unit, &
     662             :                 FMT="(T2,A,T42,F12.6,T68,F12.6,/)") &
     663        1925 :             "# Total charge", totsumorbpop(1), totsumorbpop(3)
     664             :       ELSE
     665             :          WRITE (UNIT=output_unit, &
     666             :                 FMT="(T2,A,T28,4(1X,F12.6),/)") &
     667         374 :             "# Total charge and spin", totsumorbpop(1:2), totsumorbpop(3), &
     668         748 :             totsumorbpop(1) - totsumorbpop(2)
     669             :       END IF
     670             : 
     671        2299 :       IF (output_unit > 0) CALL m_flush(output_unit)
     672             : 
     673        2299 :       CALL timestop(handle)
     674             : 
     675        2299 :    END SUBROUTINE write_orbpop
     676             : 
     677             : END MODULE population_analyses

Generated by: LCOV version 1.15