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 Routine to return block diagonal density matrix. Blocks correspond to the atomic densities
      10             : !> \par History
      11             : !>       2006.03 Moved here from qs_scf.F [Joost VandeVondele]
      12             : !>       2022.05 split from qs_initial_guess.F to break circular dependency [Harald Forbert]
      13             : ! **************************************************************************************************
      14             : MODULE qs_atomic_block
      15             :    USE atom_kind_orbitals,              ONLY: calculate_atomic_orbitals
      16             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      17             :                                               get_atomic_kind_set
      18             :    USE cp_dbcsr_api,                    ONLY: &
      19             :         dbcsr_get_info, dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, &
      20             :         dbcsr_iterator_start, dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_p_type, dbcsr_scale, &
      21             :         dbcsr_set, dbcsr_type
      22             :    USE cp_dbcsr_contrib,                ONLY: dbcsr_add_on_diag,&
      23             :                                               dbcsr_dot
      24             :    USE kinds,                           ONLY: dp
      25             :    USE message_passing,                 ONLY: mp_para_env_type
      26             :    USE qs_kind_types,                   ONLY: qs_kind_type
      27             : #include "./base/base_uses.f90"
      28             : 
      29             :    IMPLICIT NONE
      30             : 
      31             :    PRIVATE
      32             : 
      33             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_atomic_block'
      34             : 
      35             :    PUBLIC ::  calculate_atomic_block_dm
      36             : 
      37             :    TYPE atom_matrix_type
      38             :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER   :: mat => NULL()
      39             :    END TYPE atom_matrix_type
      40             : 
      41             : CONTAINS
      42             : 
      43             : ! **************************************************************************************************
      44             : !> \brief returns a block diagonal density matrix. Blocks correspond to the atomic densities.
      45             : !> \param pmatrix ...
      46             : !> \param matrix_s ...
      47             : !> \param atomic_kind_set ...
      48             : !> \param qs_kind_set ...
      49             : !> \param nspin ...
      50             : !> \param nelectron_spin ...
      51             : !> \param ounit ...
      52             : !> \param para_env ...
      53             : ! **************************************************************************************************
      54        4607 :    SUBROUTINE calculate_atomic_block_dm(pmatrix, matrix_s, atomic_kind_set, &
      55        4607 :                                         qs_kind_set, nspin, nelectron_spin, ounit, para_env)
      56             :       TYPE(dbcsr_p_type), DIMENSION(:), INTENT(INOUT)    :: pmatrix
      57             :       TYPE(dbcsr_type), INTENT(INOUT)                    :: matrix_s
      58             :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      59             :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      60             :       INTEGER, INTENT(IN)                                :: nspin
      61             :       INTEGER, DIMENSION(:), INTENT(IN)                  :: nelectron_spin
      62             :       INTEGER, INTENT(IN)                                :: ounit
      63             :       TYPE(mp_para_env_type)                             :: para_env
      64             : 
      65             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_atomic_block_dm'
      66             : 
      67             :       INTEGER                                            :: handle, icol, ikind, irow, ispin, nc, &
      68             :                                                             nkind, nocc(2)
      69        4607 :       INTEGER, ALLOCATABLE, DIMENSION(:)                 :: kind_of
      70        4607 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: nok
      71        4607 :       REAL(dp), DIMENSION(:, :), POINTER                 :: pdata
      72             :       REAL(KIND=dp)                                      :: rds, rscale, trps1
      73        4607 :       TYPE(atom_matrix_type), ALLOCATABLE, DIMENSION(:)  :: pmat
      74             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
      75             :       TYPE(dbcsr_iterator_type)                          :: iter
      76             :       TYPE(dbcsr_type), POINTER                          :: matrix_p
      77             :       TYPE(qs_kind_type), POINTER                        :: qs_kind
      78             : 
      79        4607 :       CALL timeset(routineN, handle)
      80             : 
      81        4607 :       nkind = SIZE(atomic_kind_set)
      82        4607 :       CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
      83       22127 :       ALLOCATE (pmat(nkind))
      84       13821 :       ALLOCATE (nok(2, nkind))
      85             : 
      86             :       ! precompute the atomic blocks corresponding to spherical atoms
      87       12913 :       DO ikind = 1, nkind
      88        8306 :          atomic_kind => atomic_kind_set(ikind)
      89        8306 :          qs_kind => qs_kind_set(ikind)
      90        8306 :          NULLIFY (pmat(ikind)%mat)
      91        8306 :          IF (ounit > 0) THEN
      92             :             WRITE (UNIT=ounit, FMT="(/,T2,A)") &
      93        2011 :                "Guess for atomic kind: "//TRIM(atomic_kind%name)
      94             :          END IF
      95             :          CALL calculate_atomic_orbitals(atomic_kind, qs_kind, iunit=ounit, &
      96        8306 :                                         pmat=pmat(ikind)%mat, nocc=nocc)
      97       29525 :          nok(1:2, ikind) = nocc(1:2)
      98             :       END DO
      99             : 
     100        4607 :       rscale = 1.0_dp
     101        4607 :       IF (nspin == 2) rscale = 0.5_dp
     102             : 
     103       10173 :       DO ispin = 1, nspin
     104        5566 :          IF ((ounit > 0) .AND. (nspin > 1)) THEN
     105         504 :             WRITE (UNIT=ounit, FMT="(/,T2,A,I0)") "Spin ", ispin
     106             :          END IF
     107             : 
     108        5566 :          matrix_p => pmatrix(ispin)%matrix
     109        5566 :          CALL dbcsr_set(matrix_p, 0.0_dp)
     110             : 
     111        5566 :          nocc(ispin) = 0
     112        5566 :          CALL dbcsr_iterator_start(iter, matrix_p)
     113       83489 :          DO WHILE (dbcsr_iterator_blocks_left(iter))
     114       77923 :             CALL dbcsr_iterator_next_block(iter, irow, icol, pdata)
     115       77923 :             ikind = kind_of(irow)
     116       83489 :             IF (icol .EQ. irow) THEN
     117       11784 :                IF (ispin == 1) THEN
     118             :                   pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale + &
     119     1017946 :                                 pmat(ikind)%mat(:, :, 2)*rscale
     120             :                ELSE
     121             :                   pdata(:, :) = pmat(ikind)%mat(:, :, 1)*rscale - &
     122      174908 :                                 pmat(ikind)%mat(:, :, 2)*rscale
     123             :                END IF
     124       11784 :                nocc(ispin) = nocc(ispin) + nok(ispin, ikind)
     125             :             END IF
     126             :          END DO
     127        5566 :          CALL dbcsr_iterator_stop(iter)
     128             : 
     129        5566 :          CALL dbcsr_dot(matrix_p, matrix_s, trps1)
     130        5566 :          rds = 0.0_dp
     131             :          ! could be a ghost-atoms-only simulation
     132        5566 :          IF (nelectron_spin(ispin) > 0) THEN
     133        5434 :             rds = REAL(nelectron_spin(ispin), dp)/trps1
     134             :          END IF
     135        5566 :          CALL dbcsr_scale(matrix_p, rds)
     136             : 
     137        5566 :          IF (ounit > 0) THEN
     138        1386 :             IF (nspin > 1) THEN
     139             :                WRITE (UNIT=ounit, FMT="(T2,A,I1)") &
     140         504 :                   "Re-scaling the density matrix to get the right number of electrons for spin ", ispin
     141             :             ELSE
     142             :                WRITE (UNIT=ounit, FMT="(T2,A)") &
     143         882 :                   "Re-scaling the density matrix to get the right number of electrons"
     144             :             END IF
     145        1386 :             WRITE (ounit, '(T19,A,T44,A,T67,A)') "# Electrons", "Trace(P)", "Scaling factor"
     146        1386 :             WRITE (ounit, '(T20,I10,T40,F12.3,T67,F14.3)') nelectron_spin(ispin), trps1, rds
     147             :          END IF
     148             : 
     149       21305 :          IF (nspin > 1) THEN
     150        1918 :             CALL para_env%sum(nocc)
     151        1918 :             IF (nelectron_spin(ispin) > nocc(ispin)) THEN
     152           2 :                rds = 0.99_dp
     153           2 :                CALL dbcsr_scale(matrix_p, rds)
     154           2 :                rds = (1.0_dp - rds)*nelectron_spin(ispin)
     155           2 :                CALL dbcsr_get_info(matrix_p, nfullcols_total=nc)
     156           2 :                rds = rds/REAL(nc, KIND=dp)
     157           2 :                CALL dbcsr_add_on_diag(matrix_p, rds)
     158           2 :                IF (ounit > 0) THEN
     159             :                   WRITE (UNIT=ounit, FMT="(T4,A,/,T4,A,T59,F20.12)") &
     160           1 :                      "More MOs than initial guess orbitals detected", &
     161           2 :                      "Add constant to diagonal elements ", rds
     162             :                END IF
     163             :             END IF
     164             :          END IF
     165             : 
     166             :       END DO
     167             : 
     168       12913 :       DO ikind = 1, nkind
     169       12913 :          IF (ASSOCIATED(pmat(ikind)%mat)) THEN
     170        8306 :             DEALLOCATE (pmat(ikind)%mat)
     171             :          END IF
     172             :       END DO
     173        4607 :       DEALLOCATE (pmat)
     174             : 
     175        4607 :       DEALLOCATE (kind_of, nok)
     176             : 
     177        4607 :       CALL timestop(handle)
     178             : 
     179        4607 :    END SUBROUTINE calculate_atomic_block_dm
     180             : 
     181           0 : END MODULE qs_atomic_block