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

Generated by: LCOV version 2.0-1