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