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