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
      10             : !> \par History
      11             : !>      12.2016 created [JGH]
      12             : !> \author JGH
      13             : ! **************************************************************************************************
      14             : MODULE minbas_methods
      15             :    USE cp_blacs_env,                    ONLY: cp_blacs_env_type
      16             :    USE cp_control_types,                ONLY: dft_control_type
      17             :    USE cp_dbcsr_api,                    ONLY: &
      18             :         dbcsr_create, dbcsr_distribution_type, dbcsr_filter, dbcsr_iterator_blocks_left, &
      19             :         dbcsr_iterator_next_block, dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, &
      20             :         dbcsr_multiply, dbcsr_p_type, dbcsr_release, dbcsr_reserve_diag_blocks, dbcsr_type, &
      21             :         dbcsr_type_no_symmetry
      22             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      23             :                                               copy_fm_to_dbcsr,&
      24             :                                               dbcsr_allocate_matrix_set,&
      25             :                                               dbcsr_deallocate_matrix_set
      26             :    USE cp_fm_basic_linalg,              ONLY: cp_fm_column_scale
      27             :    USE cp_fm_diag,                      ONLY: choose_eigv_solver,&
      28             :                                               cp_fm_power
      29             :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      30             :                                               cp_fm_struct_release,&
      31             :                                               cp_fm_struct_type
      32             :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      33             :                                               cp_fm_get_diag,&
      34             :                                               cp_fm_release,&
      35             :                                               cp_fm_to_fm_submat,&
      36             :                                               cp_fm_type
      37             :    USE kinds,                           ONLY: dp
      38             :    USE mao_basis,                       ONLY: mao_generate_basis
      39             :    USE message_passing,                 ONLY: mp_para_env_type
      40             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      41             :    USE particle_methods,                ONLY: get_particle_set
      42             :    USE particle_types,                  ONLY: particle_type
      43             :    USE qs_environment_types,            ONLY: get_qs_env,&
      44             :                                               qs_environment_type
      45             :    USE qs_kind_types,                   ONLY: qs_kind_type
      46             :    USE qs_ks_types,                     ONLY: get_ks_env,&
      47             :                                               qs_ks_env_type
      48             :    USE qs_mo_types,                     ONLY: get_mo_set,&
      49             :                                               mo_set_type
      50             : #include "./base/base_uses.f90"
      51             : 
      52             :    IMPLICIT NONE
      53             :    PRIVATE
      54             : 
      55             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'minbas_methods'
      56             : 
      57             :    PUBLIC ::  minbas_calculation
      58             : 
      59             : ! **************************************************************************************************
      60             : 
      61             : CONTAINS
      62             : 
      63             : ! **************************************************************************************************
      64             : !> \brief ...
      65             : !> \param qs_env ...
      66             : !> \param mos ...
      67             : !> \param quambo ...
      68             : !> \param mao ...
      69             : !> \param iounit ...
      70             : !> \param full_ortho ...
      71             : !> \param eps_filter ...
      72             : ! **************************************************************************************************
      73           0 :    SUBROUTINE minbas_calculation(qs_env, mos, quambo, mao, iounit, full_ortho, eps_filter)
      74             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      75             :       TYPE(mo_set_type), DIMENSION(:), INTENT(IN)        :: mos
      76             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: quambo
      77             :       TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
      78             :          POINTER                                         :: mao
      79             :       INTEGER, INTENT(IN), OPTIONAL                      :: iounit
      80             :       LOGICAL, INTENT(IN), OPTIONAL                      :: full_ortho
      81             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: eps_filter
      82             : 
      83             :       CHARACTER(len=*), PARAMETER :: routineN = 'minbas_calculation'
      84             : 
      85             :       INTEGER                                            :: handle, homo, i, iab, ispin, nao, natom, &
      86             :                                                             ndep, nmao, nmo, nmx, np, np1, nspin, &
      87             :                                                             nvirt, unit_nr
      88           0 :       INTEGER, DIMENSION(:), POINTER                     :: col_blk_sizes, row_blk_sizes
      89             :       LOGICAL                                            :: do_minbas, my_full_ortho
      90             :       REAL(KIND=dp)                                      :: my_eps_filter
      91           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: dval, dvalo, dvalv, eigval
      92             :       TYPE(cp_blacs_env_type), POINTER                   :: blacs_env
      93             :       TYPE(cp_fm_struct_type), POINTER                   :: fm_struct_a, fm_struct_b, fm_struct_c, &
      94             :                                                             fm_struct_d, fm_struct_e
      95             :       TYPE(cp_fm_type)                                   :: fm1, fm2, fm3, fm4, fm5, fm6, fma, fmb, &
      96             :                                                             fmwork
      97             :       TYPE(cp_fm_type), POINTER                          :: fm_mos
      98             :       TYPE(dbcsr_distribution_type), POINTER             :: dbcsr_dist
      99           0 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: mao_coef
     100           0 :       TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER       :: matrix_s
     101             :       TYPE(dbcsr_type)                                   :: smao, sortho
     102             :       TYPE(dbcsr_type), POINTER                          :: smat
     103             :       TYPE(dft_control_type), POINTER                    :: dft_control
     104             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     105           0 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     106           0 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     107             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     108             : 
     109           0 :       CALL timeset(routineN, handle)
     110             : 
     111           0 :       IF (PRESENT(iounit)) THEN
     112           0 :          unit_nr = iounit
     113             :       ELSE
     114             :          unit_nr = -1
     115             :       END IF
     116             : 
     117           0 :       IF (PRESENT(full_ortho)) THEN
     118           0 :          my_full_ortho = full_ortho
     119             :       ELSE
     120             :          my_full_ortho = .FALSE.
     121             :       END IF
     122             : 
     123           0 :       IF (PRESENT(eps_filter)) THEN
     124           0 :          my_eps_filter = eps_filter
     125             :       ELSE
     126           0 :          my_eps_filter = 1.0e-10_dp
     127             :       END IF
     128             : 
     129           0 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     130           0 :       nspin = dft_control%nspins
     131             : 
     132           0 :       CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
     133           0 :       CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
     134           0 :       CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, dbcsr_dist=dbcsr_dist)
     135           0 :       ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
     136           0 :       CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes)
     137           0 :       CALL get_particle_set(particle_set, qs_kind_set, nmao=col_blk_sizes)
     138           0 :       nmao = SUM(col_blk_sizes)
     139             :       ! check if MAOs have been specified
     140           0 :       DO iab = 1, natom
     141           0 :          IF (col_blk_sizes(iab) < 0) &
     142           0 :             CPABORT("Number of MAOs has to be specified in KIND section for all elements")
     143             :       END DO
     144           0 :       CALL get_mo_set(mo_set=mos(1), nao=nao, nmo=nmo)
     145             : 
     146           0 :       IF (unit_nr > 0) THEN
     147           0 :          WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Atomic Basis Set Functions   :', nao
     148           0 :          WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Minimal Basis Set Functions  :', nmao
     149           0 :          IF (nspin == 1) THEN
     150           0 :             WRITE (unit_nr, '(T2,A,T71,I10)') 'Total Number of Molecular Orbitals available :', nmo
     151             :          ELSE
     152           0 :             DO ispin = 1, nspin
     153           0 :                CALL get_mo_set(mo_set=mos(ispin), nmo=nmx)
     154             :                WRITE (unit_nr, '(T2,A,i2,T71,I10)') &
     155           0 :                   'Total Number of Molecular Orbitals available for Spin ', ispin, nmx
     156             :             END DO
     157             :          END IF
     158             :       END IF
     159           0 :       CPASSERT(nmao <= nao)
     160           0 :       DO ispin = 1, nspin
     161           0 :          CALL get_mo_set(mo_set=mos(ispin), nmo=nmx)
     162           0 :          IF (nmx /= nmo) EXIT
     163             :       END DO
     164           0 :       do_minbas = .TRUE.
     165           0 :       IF (nmao > nmo) THEN
     166           0 :          IF (unit_nr > 0) THEN
     167           0 :             WRITE (unit_nr, '(T2,A)') 'Localized Minimal Basis Analysis not possible'
     168             :          END IF
     169             :          do_minbas = .FALSE.
     170           0 :       ELSEIF (nmo /= nmx) THEN
     171           0 :          IF (unit_nr > 0) THEN
     172           0 :             WRITE (unit_nr, '(T2,A)') 'Different Number of Alpha and Beta MOs'
     173           0 :             WRITE (unit_nr, '(T2,A)') 'Localized Minimal Basis Analysis not possible'
     174             :          END IF
     175             :          do_minbas = .FALSE.
     176             :       ELSE
     177           0 :          IF (nao > nmo) THEN
     178           0 :             IF (unit_nr > 0) THEN
     179           0 :                WRITE (unit_nr, '(T2,A)') 'WARNING: Only a subset of MOs is available: Analysis depends on MOs'
     180             :             END IF
     181             :          END IF
     182             :       END IF
     183             : 
     184           0 :       IF (do_minbas) THEN
     185             :          ! initialize QUAMBOs
     186           0 :          NULLIFY (quambo)
     187           0 :          CALL dbcsr_allocate_matrix_set(quambo, nspin)
     188           0 :          DO ispin = 1, nspin
     189             :             ! coeficients
     190           0 :             ALLOCATE (quambo(ispin)%matrix)
     191             :             CALL dbcsr_create(matrix=quambo(ispin)%matrix, &
     192             :                               name="QUAMBO", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
     193           0 :                               row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, nze=0)
     194             :          END DO
     195             : 
     196             :          ! initialize MAOs
     197             :          ! optimize MAOs (mao_coef is allocated in the routine)
     198           0 :          CALL mao_generate_basis(qs_env, mao_coef)
     199             : 
     200             :          ! sortho (nmao x nmao)
     201             :          CALL dbcsr_create(sortho, name="SORTHO", dist=dbcsr_dist, matrix_type=dbcsr_type_no_symmetry, &
     202           0 :                            row_blk_size=col_blk_sizes, col_blk_size=col_blk_sizes, nze=0)
     203           0 :          CALL dbcsr_reserve_diag_blocks(matrix=sortho)
     204             : 
     205           0 :          DEALLOCATE (row_blk_sizes, col_blk_sizes)
     206             : 
     207             :          ! temporary FM matrices
     208           0 :          CALL get_qs_env(qs_env=qs_env, para_env=para_env, blacs_env=blacs_env)
     209           0 :          NULLIFY (fm_struct_a, fm_struct_b)
     210             :          CALL cp_fm_struct_create(fm_struct_a, nrow_global=nao, ncol_global=nmao, &
     211           0 :                                   para_env=para_env, context=blacs_env)
     212             :          CALL cp_fm_struct_create(fm_struct_b, nrow_global=nmo, ncol_global=nmao, &
     213           0 :                                   para_env=para_env, context=blacs_env)
     214           0 :          CALL cp_fm_create(fm1, fm_struct_a)
     215           0 :          CALL cp_fm_create(fm2, fm_struct_b)
     216           0 :          CALL cp_fm_create(fma, fm_struct_b)
     217           0 :          CALL cp_fm_create(fmb, fm_struct_b)
     218             : 
     219           0 :          CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
     220           0 :          smat => matrix_s(1, 1)%matrix
     221           0 :          DO ispin = 1, nspin
     222             : 
     223             :             ! SMAO = Overlap*MAOs
     224           0 :             CALL dbcsr_create(smao, name="S*MAO", template=mao_coef(1)%matrix)
     225           0 :             CALL dbcsr_multiply("N", "N", 1.0_dp, smat, mao_coef(ispin)%matrix, 0.0_dp, smao)
     226             :             ! a(nj)* = C(vn)(T) * SMAO(vj)
     227           0 :             CALL copy_dbcsr_to_fm(smao, fm1)
     228           0 :             CALL get_mo_set(mos(ispin), mo_coeff=fm_mos)
     229           0 :             CALL parallel_gemm("T", "N", nmo, nmao, nao, 1.0_dp, fm_mos, fm1, 0.0_dp, fm2)
     230           0 :             CALL dbcsr_release(smao)
     231           0 :             CALL get_mo_set(mo_set=mos(ispin), homo=homo)
     232           0 :             IF (unit_nr > 0) THEN
     233           0 :                WRITE (unit_nr, '(T2,A,T51,A,i2,T71,I10)') 'MOs in Occupied Valence Set', 'Spin ', ispin, homo
     234             :             END IF
     235           0 :             nvirt = nmo - homo
     236           0 :             NULLIFY (fm_struct_c)
     237             :             CALL cp_fm_struct_create(fm_struct_c, nrow_global=nvirt, ncol_global=nvirt, &
     238           0 :                                      para_env=para_env, context=blacs_env)
     239           0 :             CALL cp_fm_create(fm3, fm_struct_c)
     240           0 :             CALL cp_fm_create(fm4, fm_struct_c)
     241             :             ! B(vw) = a(vj)* * a(wj)*
     242             :             CALL parallel_gemm("N", "T", nvirt, nvirt, nmao, 1.0_dp, fm2, fm2, 0.0_dp, fm3, &
     243           0 :                                a_first_row=homo + 1, b_first_row=homo + 1)
     244           0 :             ALLOCATE (eigval(nvirt))
     245           0 :             CALL choose_eigv_solver(fm3, fm4, eigval)
     246             :             ! SVD(B) -> select p largest eigenvalues and vectors
     247           0 :             np = nmao - homo
     248           0 :             np1 = nvirt - np + 1
     249           0 :             IF (unit_nr > 0) THEN
     250           0 :                WRITE (unit_nr, '(T2,A,T51,A,i2,T71,I10)') 'MOs in Virtual Valence Set', 'Spin ', ispin, np
     251             :             END IF
     252             :             ! R(vw) = SUM_p T(vp)*T(wp)
     253             :             CALL parallel_gemm("N", "T", nvirt, nvirt, np, 1.0_dp, fm4, fm4, 0.0_dp, fm3, &
     254           0 :                                a_first_col=np1, b_first_col=np1)
     255             :             !
     256           0 :             ALLOCATE (dval(nmao), dvalo(nmao), dvalv(nmao))
     257           0 :             NULLIFY (fm_struct_d)
     258             :             CALL cp_fm_struct_create(fm_struct_d, nrow_global=nvirt, ncol_global=nmao, &
     259           0 :                                      para_env=para_env, context=blacs_env)
     260           0 :             CALL cp_fm_create(fm5, fm_struct_d)
     261           0 :             NULLIFY (fm_struct_e)
     262             :             CALL cp_fm_struct_create(fm_struct_e, nrow_global=nmao, ncol_global=nmao, &
     263           0 :                                      para_env=para_env, context=blacs_env)
     264           0 :             CALL cp_fm_create(fm6, fm_struct_e)
     265             :             ! D(j) = SUM_n (a(nj)*)^2 + SUM_vw R(vw) * a(vj)* * a(wj)*
     266             :             CALL parallel_gemm("N", "N", nvirt, nmao, nvirt, 1.0_dp, fm3, fm2, 0.0_dp, fm5, &
     267           0 :                                b_first_row=homo + 1)
     268             :             CALL parallel_gemm("T", "N", nmao, nmao, nvirt, 1.0_dp, fm2, fm5, 0.0_dp, fm6, &
     269           0 :                                a_first_row=homo + 1)
     270           0 :             CALL cp_fm_get_diag(fm6, dvalv(1:nmao))
     271           0 :             CALL parallel_gemm("T", "N", nmao, nmao, homo, 1.0_dp, fm2, fm2, 0.0_dp, fm6)
     272           0 :             CALL cp_fm_get_diag(fm6, dvalo(1:nmao))
     273           0 :             DO i = 1, nmao
     274           0 :                dval(i) = 1.0_dp/SQRT(dvalo(i) + dvalv(i))
     275             :             END DO
     276             :             ! scale intermediate expansion
     277           0 :             CALL cp_fm_to_fm_submat(fm2, fma, homo, nmao, 1, 1, 1, 1)
     278           0 :             CALL cp_fm_to_fm_submat(fm5, fma, nvirt, nmao, 1, 1, homo + 1, 1)
     279           0 :             CALL cp_fm_column_scale(fma, dval)
     280             :             ! Orthogonalization
     281           0 :             CALL cp_fm_create(fmwork, fm_struct_e)
     282           0 :             CALL parallel_gemm("T", "N", nmao, nmao, nmo, 1.0_dp, fma, fma, 0.0_dp, fm6)
     283           0 :             IF (my_full_ortho) THEN
     284             :                ! full orthogonalization
     285           0 :                CALL cp_fm_power(fm6, fmwork, -0.5_dp, 1.0e-12_dp, ndep)
     286           0 :                IF (ndep > 0 .AND. unit_nr > 0) THEN
     287           0 :                   WRITE (unit_nr, '(T2,A,T71,I10)') 'Warning: linear dependent basis   ', ndep
     288             :                END IF
     289           0 :                CALL parallel_gemm("N", "N", nmo, nmao, nmao, 1.0_dp, fma, fm6, 0.0_dp, fmb)
     290             :             ELSE
     291             :                ! orthogonalize on-atom blocks
     292           0 :                CALL copy_fm_to_dbcsr(fm6, sortho, keep_sparsity=.TRUE.)
     293           0 :                CALL diag_sqrt_invert(sortho)
     294           0 :                CALL copy_dbcsr_to_fm(sortho, fm6)
     295           0 :                CALL parallel_gemm("N", "N", nmo, nmao, nmao, 1.0_dp, fma, fm6, 0.0_dp, fmb)
     296             :             END IF
     297             :             ! store as QUAMBO
     298           0 :             CALL parallel_gemm("N", "N", nao, nmao, nmo, 1.0_dp, fm_mos, fmb, 0.0_dp, fm1)
     299           0 :             CALL copy_fm_to_dbcsr(fm1, quambo(ispin)%matrix)
     300           0 :             CALL dbcsr_filter(quambo(ispin)%matrix, my_eps_filter)
     301             :             !
     302           0 :             DEALLOCATE (eigval, dval, dvalo, dvalv)
     303           0 :             CALL cp_fm_release(fm3)
     304           0 :             CALL cp_fm_release(fm4)
     305           0 :             CALL cp_fm_release(fm5)
     306           0 :             CALL cp_fm_release(fm6)
     307           0 :             CALL cp_fm_release(fmwork)
     308           0 :             CALL cp_fm_struct_release(fm_struct_c)
     309           0 :             CALL cp_fm_struct_release(fm_struct_d)
     310           0 :             CALL cp_fm_struct_release(fm_struct_e)
     311             : 
     312             :          END DO
     313             : 
     314             :          ! clean up
     315           0 :          CALL cp_fm_release(fm1)
     316           0 :          CALL cp_fm_release(fm2)
     317           0 :          CALL cp_fm_release(fma)
     318           0 :          CALL cp_fm_release(fmb)
     319           0 :          CALL cp_fm_struct_release(fm_struct_a)
     320           0 :          CALL cp_fm_struct_release(fm_struct_b)
     321           0 :          CALL dbcsr_release(sortho)
     322             : 
     323             :          ! return MAOs if requested
     324           0 :          IF (PRESENT(mao)) THEN
     325           0 :             mao => mao_coef
     326             :          ELSE
     327           0 :             CALL dbcsr_deallocate_matrix_set(mao_coef)
     328             :          END IF
     329             : 
     330             :       ELSE
     331           0 :          NULLIFY (quambo)
     332             :       END IF
     333             : 
     334           0 :       CALL timestop(handle)
     335             : 
     336           0 :    END SUBROUTINE minbas_calculation
     337             : 
     338             : ! **************************************************************************************************
     339             : !> \brief ...
     340             : !> \param sortho ...
     341             : ! **************************************************************************************************
     342           0 :    SUBROUTINE diag_sqrt_invert(sortho)
     343             :       TYPE(dbcsr_type)                                   :: sortho
     344             : 
     345             :       INTEGER                                            :: i, iatom, info, jatom, lwork, n
     346           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: w, work
     347           0 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: amat, bmat
     348           0 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: sblock
     349             :       TYPE(dbcsr_iterator_type)                          :: dbcsr_iter
     350             : 
     351           0 :       CALL dbcsr_iterator_start(dbcsr_iter, sortho)
     352           0 :       DO WHILE (dbcsr_iterator_blocks_left(dbcsr_iter))
     353           0 :          CALL dbcsr_iterator_next_block(dbcsr_iter, iatom, jatom, sblock)
     354           0 :          CPASSERT(iatom == jatom)
     355           0 :          n = SIZE(sblock, 1)
     356           0 :          lwork = MAX(n*n, 100)
     357           0 :          ALLOCATE (amat(n, n), bmat(n, n), w(n), work(lwork))
     358           0 :          amat(1:n, 1:n) = sblock(1:n, 1:n)
     359           0 :          info = 0
     360           0 :          CALL dsyev("V", "U", n, amat, n, w, work, lwork, info)
     361           0 :          CPASSERT(info == 0)
     362           0 :          w(1:n) = 1._dp/SQRT(w(1:n))
     363           0 :          DO i = 1, n
     364           0 :             bmat(1:n, i) = amat(1:n, i)*w(i)
     365             :          END DO
     366           0 :          sblock(1:n, 1:n) = MATMUL(amat, TRANSPOSE(bmat))
     367           0 :          DEALLOCATE (amat, bmat, w, work)
     368             :       END DO
     369           0 :       CALL dbcsr_iterator_stop(dbcsr_iter)
     370             : 
     371           0 :    END SUBROUTINE diag_sqrt_invert
     372             : 
     373             : END MODULE minbas_methods