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             : MODULE qs_gapw_densities
      10             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      11             :                                               get_atomic_kind
      12             :    USE cp_control_types,                ONLY: dft_control_type,&
      13             :                                               gapw_control_type
      14             :    USE cp_log_handling,                 ONLY: cp_logger_get_default_io_unit
      15             :    USE kinds,                           ONLY: dp
      16             :    USE message_passing,                 ONLY: mp_para_env_type
      17             :    USE pw_env_types,                    ONLY: pw_env_get,&
      18             :                                               pw_env_type
      19             :    USE pw_pool_types,                   ONLY: pw_pool_p_type
      20             :    USE qs_charges_types,                ONLY: qs_charges_type
      21             :    USE qs_environment_types,            ONLY: get_qs_env,&
      22             :                                               qs_environment_type
      23             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      24             :                                               qs_kind_type
      25             :    USE qs_local_rho_types,              ONLY: local_rho_type
      26             :    USE qs_rho0_ggrid,                   ONLY: put_rho0_on_grid
      27             :    USE qs_rho0_methods,                 ONLY: calculate_rho0_atom
      28             :    USE qs_rho0_types,                   ONLY: rho0_atom_type,&
      29             :                                               rho0_mpole_type
      30             :    USE qs_rho_atom_methods,             ONLY: calculate_rho_atom
      31             :    USE qs_rho_atom_types,               ONLY: rho_atom_type
      32             :    USE realspace_grid_types,            ONLY: realspace_grid_desc_p_type,&
      33             :                                               realspace_grid_type
      34             : #include "./base/base_uses.f90"
      35             : 
      36             :    IMPLICIT NONE
      37             : 
      38             :    PRIVATE
      39             : 
      40             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gapw_densities'
      41             : 
      42             :    PUBLIC :: prepare_gapw_den
      43             : 
      44             : CONTAINS
      45             : 
      46             : ! **************************************************************************************************
      47             : !> \brief ...
      48             : !> \param qs_env ...
      49             : !> \param local_rho_set ...
      50             : !> \param do_rho0 ...
      51             : !> \param kind_set_external can be provided to use different projectors/grids/basis than the default
      52             : !> \param pw_env_sub ...
      53             : ! **************************************************************************************************
      54       22916 :    SUBROUTINE prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external, pw_env_sub)
      55             : 
      56             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      57             :       TYPE(local_rho_type), OPTIONAL, POINTER            :: local_rho_set
      58             :       LOGICAL, INTENT(IN), OPTIONAL                      :: do_rho0
      59             :       TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
      60             :          POINTER                                         :: kind_set_external
      61             :       TYPE(pw_env_type), OPTIONAL                        :: pw_env_sub
      62             : 
      63             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'prepare_gapw_den'
      64             : 
      65             :       INTEGER                                            :: handle, ikind, ispin, natom, nspins, &
      66             :                                                             output_unit
      67       22916 :       INTEGER, DIMENSION(:), POINTER                     :: atom_list
      68             :       LOGICAL                                            :: extern, my_do_rho0, paw_atom
      69             :       REAL(dp)                                           :: rho0_h_tot, tot_rs_int
      70       22916 :       REAL(dp), DIMENSION(:), POINTER                    :: rho1_h_tot, rho1_s_tot
      71       22916 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
      72             :       TYPE(dft_control_type), POINTER                    :: dft_control
      73             :       TYPE(gapw_control_type), POINTER                   :: gapw_control
      74             :       TYPE(mp_para_env_type), POINTER                    :: para_env
      75       22916 :       TYPE(pw_pool_p_type), DIMENSION(:), POINTER        :: my_pools
      76             :       TYPE(qs_charges_type), POINTER                     :: qs_charges
      77       22916 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: my_kind_set
      78             :       TYPE(realspace_grid_desc_p_type), DIMENSION(:), &
      79       22916 :          POINTER                                         :: my_rs_descs
      80       22916 :       TYPE(realspace_grid_type), DIMENSION(:), POINTER   :: my_rs_grids
      81       22916 :       TYPE(rho0_atom_type), DIMENSION(:), POINTER        :: rho0_atom_set
      82             :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
      83       22916 :       TYPE(rho_atom_type), DIMENSION(:), POINTER         :: rho_atom_set
      84             : 
      85       22916 :       CALL timeset(routineN, handle)
      86             : 
      87       22916 :       NULLIFY (atomic_kind_set)
      88       22916 :       NULLIFY (my_kind_set)
      89       22916 :       NULLIFY (dft_control)
      90       22916 :       NULLIFY (gapw_control)
      91       22916 :       NULLIFY (para_env)
      92       22916 :       NULLIFY (atom_list)
      93       22916 :       NULLIFY (rho0_mpole)
      94       22916 :       NULLIFY (qs_charges)
      95       22916 :       NULLIFY (rho1_h_tot, rho1_s_tot)
      96       22916 :       NULLIFY (rho_atom_set)
      97       22916 :       NULLIFY (rho0_atom_set)
      98             : 
      99       22916 :       my_do_rho0 = .TRUE.
     100       22916 :       IF (PRESENT(do_rho0)) my_do_rho0 = do_rho0
     101             : 
     102       22916 :       output_unit = cp_logger_get_default_io_unit()
     103             : 
     104             :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
     105             :                       para_env=para_env, &
     106             :                       qs_charges=qs_charges, &
     107             :                       qs_kind_set=my_kind_set, &
     108             :                       atomic_kind_set=atomic_kind_set, &
     109             :                       rho0_mpole=rho0_mpole, &
     110             :                       rho_atom_set=rho_atom_set, &
     111       22916 :                       rho0_atom_set=rho0_atom_set)
     112             : 
     113       22916 :       gapw_control => dft_control%qs_control%gapw_control
     114             : 
     115             :       ! If TDDFPT%MGRID is defined, overwrite QS grid info accordingly
     116       22916 :       IF (PRESENT(local_rho_set)) THEN
     117        6190 :          rho_atom_set => local_rho_set%rho_atom_set
     118        6190 :          IF (my_do_rho0) THEN
     119        2196 :             rho0_mpole => local_rho_set%rho0_mpole
     120        2196 :             rho0_atom_set => local_rho_set%rho0_atom_set
     121             :          END IF
     122             :       END IF
     123             : 
     124       22916 :       extern = .FALSE.
     125       22916 :       IF (PRESENT(kind_set_external)) THEN
     126        2844 :          CPASSERT(ASSOCIATED(kind_set_external))
     127        2844 :          my_kind_set => kind_set_external
     128        2844 :          extern = .TRUE.
     129             :       END IF
     130             : 
     131       22916 :       nspins = dft_control%nspins
     132             : 
     133       22916 :       rho0_h_tot = 0.0_dp
     134      114580 :       ALLOCATE (rho1_h_tot(1:nspins), rho1_s_tot(1:nspins))
     135       50022 :       rho1_h_tot = 0.0_dp
     136       50022 :       rho1_s_tot = 0.0_dp
     137             : 
     138       70724 :       DO ikind = 1, SIZE(atomic_kind_set)
     139       47808 :          CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
     140       47808 :          CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
     141             : 
     142             :          !Calculate rho1_h and rho1_s on the radial grids centered on the atomic position
     143       47808 :          IF (paw_atom) THEN
     144             :             CALL calculate_rho_atom(para_env, rho_atom_set, my_kind_set(ikind), &
     145       43670 :                                     atom_list, natom, nspins, rho1_h_tot, rho1_s_tot)
     146             :          END IF
     147             : 
     148             :          !Calculate rho0_h and rho0_s on the radial grids centered on the atomic position
     149       47808 :          IF (my_do_rho0) &
     150             :             CALL calculate_rho0_atom(gapw_control, rho_atom_set, rho0_atom_set, rho0_mpole, &
     151      103294 :                                      atom_list, natom, ikind, my_kind_set(ikind), rho0_h_tot)
     152             : 
     153             :       END DO
     154             : 
     155             :       !Do not mess with charges if using a non-default kind_set
     156       22916 :       IF (.NOT. extern) THEN
     157       67256 :          CALL para_env%sum(rho1_h_tot)
     158       67256 :          CALL para_env%sum(rho1_s_tot)
     159       43664 :          DO ispin = 1, nspins
     160       23592 :             qs_charges%total_rho1_hard(ispin) = -rho1_h_tot(ispin)
     161       43664 :             qs_charges%total_rho1_soft(ispin) = -rho1_s_tot(ispin)
     162             :          END DO
     163             : 
     164       20072 :          IF (my_do_rho0) THEN
     165       16010 :             rho0_mpole%total_rho0_h = -rho0_h_tot
     166             : 
     167             :             ! When MGRID is defined within TDDFPT
     168       16010 :             IF (PRESENT(pw_env_sub)) THEN
     169             :                ! Find pool
     170         972 :                NULLIFY (my_pools, my_rs_grids, my_rs_descs)
     171             :                CALL pw_env_get(pw_env=pw_env_sub, rs_grids=my_rs_grids, &
     172         972 :                                rs_descs=my_rs_descs, pw_pools=my_pools)
     173             :                ! Put the rho0_soft on the global grid
     174             :                CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int, my_pools=my_pools, &
     175         972 :                                      my_rs_grids=my_rs_grids, my_rs_descs=my_rs_descs)
     176             :             ELSE
     177             :                ! Put the rho0_soft on the global grid
     178       15038 :                CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int)
     179             :             END IF
     180             : 
     181       16010 :             IF (ABS(rho0_h_tot) .GE. 1.0E-5_dp) THEN
     182       15084 :                IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_h_tot)) .GT. 1.0E-3_dp) THEN
     183        1778 :                   IF (output_unit > 0) THEN
     184         889 :                      WRITE (output_unit, '(/,72("*"))')
     185             :                      WRITE (output_unit, '(T2,A,T66,1E20.8)') &
     186         889 :                         "WARNING: rho0 calculated on the local grid is  :", -rho0_h_tot, &
     187        1778 :                         "         rho0 calculated on the global grid is :", tot_rs_int
     188             :                      WRITE (output_unit, '(T2,A)') &
     189         889 :                         "         bad integration"
     190         889 :                      WRITE (output_unit, '(72("*"),/)')
     191             :                   END IF
     192             :                END IF
     193             :             END IF
     194       16010 :             qs_charges%total_rho0_soft_rspace = tot_rs_int
     195       16010 :             qs_charges%total_rho0_hard_lebedev = rho0_h_tot
     196             :          ELSE
     197        4062 :             qs_charges%total_rho0_hard_lebedev = 0.0_dp
     198             :          END IF
     199             :       END IF
     200             : 
     201       22916 :       DEALLOCATE (rho1_h_tot, rho1_s_tot)
     202             : 
     203       22916 :       CALL timestop(handle)
     204             : 
     205       22916 :    END SUBROUTINE prepare_gapw_den
     206             : 
     207             : END MODULE qs_gapw_densities