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             : MODULE qs_rho0_types
       9             : 
      10             :    USE cp_units,                        ONLY: cp_unit_from_cp2k
      11             :    USE kinds,                           ONLY: dp
      12             :    USE mathconstants,                   ONLY: fourpi,&
      13             :                                               pi,&
      14             :                                               rootpi
      15             :    USE memory_utilities,                ONLY: reallocate
      16             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      17             :                                               pw_r3d_rs_type
      18             :    USE qs_grid_atom,                    ONLY: grid_atom_type
      19             :    USE qs_rho_atom_types,               ONLY: rho_atom_coeff
      20             :    USE whittaker,                       ONLY: whittaker_c0a,&
      21             :                                               whittaker_ci
      22             : #include "./base/base_uses.f90"
      23             : 
      24             :    IMPLICIT NONE
      25             : 
      26             :    PRIVATE
      27             : 
      28             : ! *** Global parameters (only in this module)
      29             : 
      30             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_rho0_types'
      31             : 
      32             : ! *** Define multipole type ***
      33             : 
      34             : ! **************************************************************************************************
      35             :    TYPE mpole_rho_atom
      36             :       REAL(dp), DIMENSION(:), POINTER             ::  Qlm_h => NULL(), &
      37             :                                                      Qlm_s => NULL(), &
      38             :                                                      Qlm_tot => NULL(), &
      39             :                                                      Qlm_car => NULL()
      40             :       REAL(dp)                                    ::  Qlm_z = -1.0_dp
      41             :       REAL(dp), DIMENSION(2)                      ::  Q0 = -1.0_dp
      42             :    END TYPE mpole_rho_atom
      43             : 
      44             : ! **************************************************************************************************
      45             :    TYPE mpole_gau_overlap
      46             :       REAL(dp), DIMENSION(:, :, :), POINTER         :: Qlm_gg => NULL()
      47             :       REAL(dp), DIMENSION(:, :), POINTER           :: g0_h => NULL(), vg0_h => NULL()
      48             :       REAL(dp)                                    :: rpgf0_h = -1.0_dp, rpgf0_s = -1.0_dp
      49             :    END TYPE mpole_gau_overlap
      50             : 
      51             : ! **************************************************************************************************
      52             :    TYPE rho0_mpole_type
      53             :       TYPE(mpole_rho_atom), DIMENSION(:), POINTER  :: mp_rho => NULL()
      54             :       TYPE(mpole_gau_overlap), DIMENSION(:), &
      55             :          POINTER   :: mp_gau => NULL()
      56             :       REAL(dp)                                    :: zet0_h = -1.0_dp, &
      57             :                                                      total_rho0_h = -1.0_dp
      58             :       REAL(dp)                                    :: max_rpgf0_s = -1.0_dp
      59             :       REAL(dp), DIMENSION(:), POINTER             :: norm_g0l_h => NULL()
      60             :       INTEGER, DIMENSION(:), POINTER             :: lmax0_kind => NULL()
      61             :       INTEGER                                     :: lmax_0 = -1, igrid_zet0_s = -1
      62             :       TYPE(pw_r3d_rs_type), POINTER                    :: rho0_s_rs => NULL()
      63             :       TYPE(pw_c1d_gs_type), POINTER ::              rho0_s_gs => NULL()
      64             :    END TYPE rho0_mpole_type
      65             : 
      66             : ! **************************************************************************************************
      67             :    TYPE rho0_atom_type
      68             :       TYPE(rho_atom_coeff), POINTER               :: rho0_rad_h => NULL(), &
      69             :                                                      vrho0_rad_h => NULL()
      70             :    END TYPE rho0_atom_type
      71             : 
      72             : ! Public Types
      73             : 
      74             :    PUBLIC :: mpole_rho_atom, mpole_gau_overlap, &
      75             :              rho0_atom_type, rho0_mpole_type
      76             : 
      77             : ! Public Subroutine
      78             : 
      79             :    PUBLIC :: allocate_multipoles, allocate_rho0_mpole, &
      80             :              allocate_rho0_atom, allocate_rho0_atom_rad, &
      81             :              deallocate_rho0_atom, deallocate_rho0_mpole, &
      82             :              calculate_g0, get_rho0_mpole, initialize_mpole_rho, &
      83             :              write_rho0_info
      84             : 
      85             : CONTAINS
      86             : 
      87             : ! **************************************************************************************************
      88             : !> \brief ...
      89             : !> \param mp_rho ...
      90             : !> \param natom ...
      91             : !> \param mp_gau ...
      92             : !> \param nkind ...
      93             : ! **************************************************************************************************
      94        1790 :    SUBROUTINE allocate_multipoles(mp_rho, natom, mp_gau, nkind)
      95             : 
      96             :       TYPE(mpole_rho_atom), DIMENSION(:), POINTER        :: mp_rho
      97             :       INTEGER, INTENT(IN)                                :: natom
      98             :       TYPE(mpole_gau_overlap), DIMENSION(:), POINTER     :: mp_gau
      99             :       INTEGER, INTENT(IN)                                :: nkind
     100             : 
     101             :       INTEGER                                            :: iat, ikind
     102             : 
     103        1790 :       IF (ASSOCIATED(mp_rho)) THEN
     104           0 :          CALL deallocate_mpole_rho(mp_rho)
     105             :       END IF
     106             : 
     107       15076 :       ALLOCATE (mp_rho(natom))
     108             : 
     109        7916 :       DO iat = 1, natom
     110        6126 :          NULLIFY (mp_rho(iat)%Qlm_h)
     111        6126 :          NULLIFY (mp_rho(iat)%Qlm_s)
     112        6126 :          NULLIFY (mp_rho(iat)%Qlm_tot)
     113        7916 :          NULLIFY (mp_rho(iat)%Qlm_car)
     114             :       END DO
     115             : 
     116        1790 :       IF (ASSOCIATED(mp_gau)) THEN
     117           0 :          CALL deallocate_mpole_gau(mp_gau)
     118             :       END IF
     119             : 
     120        9112 :       ALLOCATE (mp_gau(nkind))
     121             : 
     122        5532 :       DO ikind = 1, nkind
     123        3742 :          NULLIFY (mp_gau(ikind)%Qlm_gg)
     124        3742 :          NULLIFY (mp_gau(ikind)%g0_h)
     125        3742 :          NULLIFY (mp_gau(ikind)%vg0_h)
     126        3742 :          mp_gau(ikind)%rpgf0_h = 0.0_dp
     127        5532 :          mp_gau(ikind)%rpgf0_s = 0.0_dp
     128             :       END DO
     129             : 
     130        1790 :    END SUBROUTINE allocate_multipoles
     131             : 
     132             : ! **************************************************************************************************
     133             : !> \brief ...
     134             : !> \param rho0_set ...
     135             : !> \param natom ...
     136             : ! **************************************************************************************************
     137        1790 :    SUBROUTINE allocate_rho0_atom(rho0_set, natom)
     138             : 
     139             :       TYPE(rho0_atom_type), DIMENSION(:), POINTER        :: rho0_set
     140             :       INTEGER, INTENT(IN)                                :: natom
     141             : 
     142             :       INTEGER                                            :: iat
     143             : 
     144        1790 :       IF (ASSOCIATED(rho0_set)) THEN
     145           0 :          CALL deallocate_rho0_atom(rho0_set)
     146             :       END IF
     147             : 
     148       11496 :       ALLOCATE (rho0_set(natom))
     149             : 
     150        7916 :       DO iat = 1, natom
     151        6126 :          NULLIFY (rho0_set(iat)%rho0_rad_h)
     152        7916 :          NULLIFY (rho0_set(iat)%vrho0_rad_h)
     153             :       END DO
     154             : 
     155        1790 :    END SUBROUTINE allocate_rho0_atom
     156             : 
     157             : ! **************************************************************************************************
     158             : !> \brief ...
     159             : !> \param rho0_atom ...
     160             : !> \param nr ...
     161             : !> \param nchannels ...
     162             : ! **************************************************************************************************
     163        6126 :    SUBROUTINE allocate_rho0_atom_rad(rho0_atom, nr, nchannels)
     164             : 
     165             :       TYPE(rho0_atom_type), INTENT(OUT)                  :: rho0_atom
     166             :       INTEGER, INTENT(IN)                                :: nr, nchannels
     167             : 
     168        6126 :       ALLOCATE (rho0_atom%rho0_rad_h)
     169             : 
     170             :       NULLIFY (rho0_atom%rho0_rad_h%r_coef)
     171       24504 :       ALLOCATE (rho0_atom%rho0_rad_h%r_coef(1:nr, 1:nchannels))
     172     2393000 :       rho0_atom%rho0_rad_h%r_coef = 0.0_dp
     173             : 
     174        6126 :       ALLOCATE (rho0_atom%vrho0_rad_h)
     175             : 
     176             :       NULLIFY (rho0_atom%vrho0_rad_h%r_coef)
     177       24504 :       ALLOCATE (rho0_atom%vrho0_rad_h%r_coef(1:nr, 1:nchannels))
     178     2393000 :       rho0_atom%vrho0_rad_h%r_coef = 0.0_dp
     179             : 
     180        6126 :    END SUBROUTINE allocate_rho0_atom_rad
     181             : 
     182             : ! **************************************************************************************************
     183             : !> \brief ...
     184             : !> \param rho0 ...
     185             : ! **************************************************************************************************
     186        1790 :    SUBROUTINE allocate_rho0_mpole(rho0)
     187             : 
     188             :       TYPE(rho0_mpole_type), POINTER                     :: rho0
     189             : 
     190        1790 :       IF (ASSOCIATED(rho0)) THEN
     191           0 :          CALL deallocate_rho0_mpole(rho0)
     192             :       END IF
     193             : 
     194        1790 :       ALLOCATE (rho0)
     195             : 
     196             :       NULLIFY (rho0%mp_rho)
     197             :       NULLIFY (rho0%mp_gau)
     198             :       NULLIFY (rho0%norm_g0l_h)
     199             :       NULLIFY (rho0%lmax0_kind)
     200             :       NULLIFY (rho0%rho0_s_rs)
     201             :       NULLIFY (rho0%rho0_s_gs)
     202             : 
     203        1790 :    END SUBROUTINE allocate_rho0_mpole
     204             : 
     205             : ! **************************************************************************************************
     206             : !> \brief ...
     207             : !> \param rho0_mpole ...
     208             : !> \param grid_atom ...
     209             : !> \param ik ...
     210             : ! **************************************************************************************************
     211        3742 :    SUBROUTINE calculate_g0(rho0_mpole, grid_atom, ik)
     212             : 
     213             :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
     214             :       TYPE(grid_atom_type), POINTER                      :: grid_atom
     215             :       INTEGER, INTENT(IN)                                :: ik
     216             : 
     217             :       INTEGER                                            :: ir, l, lmax, nr
     218             :       REAL(dp)                                           :: c1, prefactor, root_z_h, z_h
     219        3742 :       REAL(dp), ALLOCATABLE, DIMENSION(:)                :: erf_z_h, gexp, gh_tmp, int1, int2
     220             : 
     221        3742 :       nr = grid_atom%nr
     222        3742 :       lmax = rho0_mpole%lmax0_kind(ik)
     223        3742 :       z_h = rho0_mpole%zet0_h
     224        3742 :       root_z_h = SQRT(z_h)
     225             : 
     226             : !   Allocate g0
     227        3742 :       CALL reallocate(rho0_mpole%mp_gau(ik)%g0_h, 1, nr, 0, lmax)
     228        3742 :       CALL reallocate(rho0_mpole%mp_gau(ik)%vg0_h, 1, nr, 0, lmax)
     229             : 
     230       26194 :       ALLOCATE (gexp(nr), gh_tmp(nr), erf_z_h(nr), int1(nr), int2(nr))
     231             : 
     232      195762 :       gh_tmp(1:nr) = EXP(-z_h*grid_atom%rad2(1:nr))
     233             : 
     234      195762 :       DO ir = 1, nr
     235      195762 :          erf_z_h(ir) = erf(grid_atom%rad(ir)*root_z_h)
     236             :       END DO
     237             : 
     238      195762 :       DO ir = 1, nr
     239      195762 :          IF (gh_tmp(ir) < 1.0E-30_dp) gh_tmp(ir) = 0.0_dp
     240             :       END DO
     241             : 
     242      195762 :       gexp(1:nr) = gh_tmp(1:nr)
     243             :       rho0_mpole%mp_gau(ik)%g0_h(1:nr, 0) = gh_tmp(1:nr)* &
     244      195762 :                                             rho0_mpole%norm_g0l_h(0)
     245        3742 :       CALL whittaker_c0a(int1, grid_atom%rad, gh_tmp, erf_z_h, z_h, 0, 0, nr)
     246        3742 :       CALL whittaker_ci(int2, grid_atom%rad, gh_tmp, z_h, 0, nr)
     247             : 
     248        3742 :       prefactor = fourpi*rho0_mpole%norm_g0l_h(0)
     249             : 
     250        3742 :       c1 = SQRT(pi*pi*pi/(z_h*z_h*z_h))*rho0_mpole%norm_g0l_h(0)
     251             : 
     252      195762 :       DO ir = 1, nr
     253      195762 :          rho0_mpole%mp_gau(ik)%vg0_h(ir, 0) = c1*erf_z_h(ir)*grid_atom%oorad2l(ir, 1)
     254             :       END DO
     255             : 
     256       10226 :       DO l = 1, lmax
     257      339644 :          gh_tmp(1:nr) = gh_tmp(1:nr)*grid_atom%rad(1:nr)
     258             :          rho0_mpole%mp_gau(ik)%g0_h(1:nr, l) = gh_tmp(1:nr)* &
     259      339644 :                                                rho0_mpole%norm_g0l_h(l)
     260             : 
     261        6484 :          prefactor = fourpi/(2.0_dp*l + 1.0_dp)*rho0_mpole%norm_g0l_h(l)
     262        6484 :          CALL whittaker_c0a(int1, grid_atom%rad, gexp, erf_z_h, z_h, l, l, nr)
     263      343386 :          DO ir = 1, nr
     264             :             rho0_mpole%mp_gau(ik)%vg0_h(ir, l) = prefactor*(int1(ir) + &
     265      339644 :                                                             int2(ir)*grid_atom%rad2l(ir, l))
     266             :          END DO
     267             : 
     268             :       END DO ! l
     269             : 
     270        3742 :       DEALLOCATE (gexp, erf_z_h, gh_tmp, int1, int2)
     271        3742 :    END SUBROUTINE calculate_g0
     272             : 
     273             : ! **************************************************************************************************
     274             : !> \brief ...
     275             : !> \param mp_gau ...
     276             : ! **************************************************************************************************
     277        1790 :    SUBROUTINE deallocate_mpole_gau(mp_gau)
     278             : 
     279             :       TYPE(mpole_gau_overlap), DIMENSION(:), POINTER     :: mp_gau
     280             : 
     281             :       INTEGER                                            :: ikind, nkind
     282             : 
     283        1790 :       nkind = SIZE(mp_gau)
     284             : 
     285        5532 :       DO ikind = 1, nkind
     286             : 
     287        3742 :          IF (ASSOCIATED(mp_gau(ikind)%Qlm_gg)) THEN
     288        3400 :             DEALLOCATE (mp_gau(ikind)%Qlm_gg)
     289             :          END IF
     290             : 
     291        3742 :          DEALLOCATE (mp_gau(ikind)%g0_h)
     292             : 
     293        5532 :          DEALLOCATE (mp_gau(ikind)%vg0_h)
     294             :       END DO
     295             : 
     296        1790 :       DEALLOCATE (mp_gau)
     297             : 
     298        1790 :    END SUBROUTINE deallocate_mpole_gau
     299             : 
     300             : ! **************************************************************************************************
     301             : !> \brief ...
     302             : !> \param mp_rho ...
     303             : ! **************************************************************************************************
     304        1790 :    SUBROUTINE deallocate_mpole_rho(mp_rho)
     305             : 
     306             :       TYPE(mpole_rho_atom), DIMENSION(:), POINTER        :: mp_rho
     307             : 
     308             :       INTEGER                                            :: iat, natom
     309             : 
     310        1790 :       natom = SIZE(mp_rho)
     311             : 
     312        7916 :       DO iat = 1, natom
     313        6126 :          DEALLOCATE (mp_rho(iat)%Qlm_h)
     314        6126 :          DEALLOCATE (mp_rho(iat)%Qlm_s)
     315        6126 :          DEALLOCATE (mp_rho(iat)%Qlm_tot)
     316        7916 :          DEALLOCATE (mp_rho(iat)%Qlm_car)
     317             :       END DO
     318             : 
     319        1790 :       DEALLOCATE (mp_rho)
     320             : 
     321        1790 :    END SUBROUTINE deallocate_mpole_rho
     322             : 
     323             : ! **************************************************************************************************
     324             : !> \brief ...
     325             : !> \param rho0_atom_set ...
     326             : ! **************************************************************************************************
     327        1790 :    SUBROUTINE deallocate_rho0_atom(rho0_atom_set)
     328             : 
     329             :       TYPE(rho0_atom_type), DIMENSION(:), POINTER        :: rho0_atom_set
     330             : 
     331             :       INTEGER                                            :: iat, natom
     332             : 
     333        1790 :       IF (ASSOCIATED(rho0_atom_set)) THEN
     334             : 
     335        1790 :          natom = SIZE(rho0_atom_set)
     336             : 
     337        7916 :          DO iat = 1, natom
     338        6126 :             IF (ASSOCIATED(rho0_atom_set(iat)%rho0_rad_h)) THEN
     339        6126 :                DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h%r_coef)
     340        6126 :                DEALLOCATE (rho0_atom_set(iat)%rho0_rad_h)
     341             :             END IF
     342        7916 :             IF (ASSOCIATED(rho0_atom_set(iat)%vrho0_rad_h)) THEN
     343        6126 :                DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h%r_coef)
     344        6126 :                DEALLOCATE (rho0_atom_set(iat)%vrho0_rad_h)
     345             :             END IF
     346             :          END DO
     347             : 
     348        1790 :          DEALLOCATE (rho0_atom_set)
     349             :       ELSE
     350             :          CALL cp_abort(__LOCATION__, &
     351             :                        "The pointer rho0_atom_set is not associated and "// &
     352           0 :                        "cannot be deallocated")
     353             :       END IF
     354             : 
     355        1790 :    END SUBROUTINE deallocate_rho0_atom
     356             : ! **************************************************************************************************
     357             : !> \brief ...
     358             : !> \param rho0 ...
     359             : ! **************************************************************************************************
     360        1790 :    SUBROUTINE deallocate_rho0_mpole(rho0)
     361             : 
     362             :       TYPE(rho0_mpole_type), POINTER                     :: rho0
     363             : 
     364        1790 :       IF (ASSOCIATED(rho0)) THEN
     365             : 
     366        1790 :          IF (ASSOCIATED(rho0%mp_gau)) CALL deallocate_mpole_gau(rho0%mp_gau)
     367             : 
     368        1790 :          IF (ASSOCIATED(rho0%mp_rho)) CALL deallocate_mpole_rho(rho0%mp_rho)
     369             : 
     370        1790 :          IF (ASSOCIATED(rho0%lmax0_kind)) THEN
     371        1790 :             DEALLOCATE (rho0%lmax0_kind)
     372             :          END IF
     373             : 
     374        1790 :          IF (ASSOCIATED(rho0%norm_g0l_h)) THEN
     375        1790 :             DEALLOCATE (rho0%norm_g0l_h)
     376             :          END IF
     377             : 
     378        1790 :          IF (ASSOCIATED(rho0%rho0_s_rs)) THEN
     379        1790 :             CALL rho0%rho0_s_rs%release()
     380        1790 :             DEALLOCATE (rho0%rho0_s_rs)
     381             :          END IF
     382             : 
     383        1790 :          IF (ASSOCIATED(rho0%rho0_s_gs)) THEN
     384        1790 :             CALL rho0%rho0_s_gs%release()
     385        1790 :             DEALLOCATE (rho0%rho0_s_gs)
     386             : 
     387             :          END IF
     388        1790 :          DEALLOCATE (rho0)
     389             :       ELSE
     390             :          CALL cp_abort(__LOCATION__, &
     391             :                        "The pointer rho0 is not associated and "// &
     392           0 :                        "cannot be deallocated")
     393             :       END IF
     394             : 
     395        1790 :    END SUBROUTINE deallocate_rho0_mpole
     396             : 
     397             : ! **************************************************************************************************
     398             : !> \brief ...
     399             : !> \param rho0_mpole ...
     400             : !> \param g0_h ...
     401             : !> \param vg0_h ...
     402             : !> \param iat ...
     403             : !> \param ikind ...
     404             : !> \param lmax_0 ...
     405             : !> \param l0_ikind ...
     406             : !> \param mp_gau_ikind ...
     407             : !> \param mp_rho ...
     408             : !> \param norm_g0l_h ...
     409             : !> \param Qlm_gg ...
     410             : !> \param Qlm_car ...
     411             : !> \param Qlm_tot ...
     412             : !> \param zet0_h ...
     413             : !> \param igrid_zet0_s ...
     414             : !> \param rpgf0_h ...
     415             : !> \param rpgf0_s ...
     416             : !> \param max_rpgf0_s ...
     417             : !> \param rho0_s_rs ...
     418             : !> \param rho0_s_gs ...
     419             : ! **************************************************************************************************
     420      221865 :    SUBROUTINE get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, &
     421             :                              mp_gau_ikind, mp_rho, norm_g0l_h, &
     422             :                              Qlm_gg, Qlm_car, Qlm_tot, &
     423             :                              zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, &
     424             :                              max_rpgf0_s, rho0_s_rs, rho0_s_gs)
     425             : 
     426             :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
     427             :       REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER       :: g0_h, vg0_h
     428             :       INTEGER, INTENT(IN), OPTIONAL                      :: iat, ikind
     429             :       INTEGER, INTENT(OUT), OPTIONAL                     :: lmax_0, l0_ikind
     430             :       TYPE(mpole_gau_overlap), OPTIONAL, POINTER         :: mp_gau_ikind
     431             :       TYPE(mpole_rho_atom), DIMENSION(:), OPTIONAL, &
     432             :          POINTER                                         :: mp_rho
     433             :       REAL(dp), DIMENSION(:), OPTIONAL, POINTER          :: norm_g0l_h
     434             :       REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER    :: Qlm_gg
     435             :       REAL(dp), DIMENSION(:), OPTIONAL, POINTER          :: Qlm_car, Qlm_tot
     436             :       REAL(dp), INTENT(OUT), OPTIONAL                    :: zet0_h
     437             :       INTEGER, INTENT(OUT), OPTIONAL                     :: igrid_zet0_s
     438             :       REAL(dp), INTENT(OUT), OPTIONAL                    :: rpgf0_h, rpgf0_s, max_rpgf0_s
     439             :       TYPE(pw_r3d_rs_type), OPTIONAL, POINTER            :: rho0_s_rs
     440             :       TYPE(pw_c1d_gs_type), OPTIONAL, POINTER            :: rho0_s_gs
     441             : 
     442      221865 :       IF (ASSOCIATED(rho0_mpole)) THEN
     443             : 
     444      221865 :          IF (PRESENT(lmax_0)) lmax_0 = rho0_mpole%lmax_0
     445      221865 :          IF (PRESENT(mp_rho)) mp_rho => rho0_mpole%mp_rho
     446      221865 :          IF (PRESENT(norm_g0l_h)) norm_g0l_h => rho0_mpole%norm_g0l_h
     447      221865 :          IF (PRESENT(zet0_h)) zet0_h = rho0_mpole%zet0_h
     448      221865 :          IF (PRESENT(igrid_zet0_s)) igrid_zet0_s = rho0_mpole%igrid_zet0_s
     449      221865 :          IF (PRESENT(max_rpgf0_s)) max_rpgf0_s = rho0_mpole%max_rpgf0_s
     450      221865 :          IF (PRESENT(rho0_s_rs)) rho0_s_rs => rho0_mpole%rho0_s_rs
     451      221865 :          IF (PRESENT(rho0_s_gs)) rho0_s_gs => rho0_mpole%rho0_s_gs
     452             : 
     453      221865 :          IF (PRESENT(ikind)) THEN
     454      125952 :             IF (PRESENT(l0_ikind)) l0_ikind = rho0_mpole%lmax0_kind(ikind)
     455      125952 :             IF (PRESENT(mp_gau_ikind)) mp_gau_ikind => rho0_mpole%mp_gau(ikind)
     456      125952 :             IF (PRESENT(g0_h)) g0_h => rho0_mpole%mp_gau(ikind)%g0_h
     457      125952 :             IF (PRESENT(vg0_h)) vg0_h => rho0_mpole%mp_gau(ikind)%vg0_h
     458      125952 :             IF (PRESENT(Qlm_gg)) Qlm_gg => rho0_mpole%mp_gau(ikind)%Qlm_gg
     459      125952 :             IF (PRESENT(rpgf0_h)) rpgf0_h = rho0_mpole%mp_gau(ikind)%rpgf0_h
     460      125952 :             IF (PRESENT(rpgf0_s)) rpgf0_s = rho0_mpole%mp_gau(ikind)%rpgf0_s
     461             :          END IF
     462      221865 :          IF (PRESENT(iat)) THEN
     463       47799 :             IF (PRESENT(Qlm_car)) Qlm_car => rho0_mpole%mp_rho(iat)%Qlm_car
     464       47799 :             IF (PRESENT(Qlm_tot)) Qlm_tot => rho0_mpole%mp_rho(iat)%Qlm_tot
     465             :          END IF
     466             : 
     467             :       ELSE
     468           0 :          CPABORT("The pointer rho0_mpole is not associated")
     469             :       END IF
     470             : 
     471      221865 :    END SUBROUTINE get_rho0_mpole
     472             : 
     473             : ! **************************************************************************************************
     474             : !> \brief ...
     475             : !> \param mp_rho ...
     476             : !> \param nchan_s ...
     477             : !> \param nchan_c ...
     478             : !> \param zeff ...
     479             : ! **************************************************************************************************
     480        6126 :    SUBROUTINE initialize_mpole_rho(mp_rho, nchan_s, nchan_c, zeff)
     481             : 
     482             :       TYPE(mpole_rho_atom)                               :: mp_rho
     483             :       INTEGER, INTENT(IN)                                :: nchan_s, nchan_c
     484             :       REAL(KIND=dp), INTENT(IN)                          :: zeff
     485             : 
     486        6126 :       CALL reallocate(mp_rho%Qlm_h, 1, nchan_s)
     487        6126 :       CALL reallocate(mp_rho%Qlm_s, 1, nchan_s)
     488        6126 :       CALL reallocate(mp_rho%Qlm_tot, 1, nchan_s)
     489        6126 :       CALL reallocate(mp_rho%Qlm_car, 1, nchan_c)
     490             : 
     491       52060 :       mp_rho%Qlm_h = 0.0_dp
     492       52060 :       mp_rho%Qlm_s = 0.0_dp
     493       52060 :       mp_rho%Qlm_tot = 0.0_dp
     494       57904 :       mp_rho%Qlm_car = 0.0_dp
     495        6126 :       mp_rho%Qlm_z = -2.0_dp*rootpi*Zeff
     496       18378 :       mp_rho%Q0 = 0.0_dp
     497             : 
     498        6126 :    END SUBROUTINE initialize_mpole_rho
     499             : 
     500             : ! **************************************************************************************************
     501             : !> \brief ...
     502             : !> \param rho0_mpole ...
     503             : !> \param unit_str ...
     504             : !> \param output_unit ...
     505             : ! **************************************************************************************************
     506           1 :    SUBROUTINE write_rho0_info(rho0_mpole, unit_str, output_unit)
     507             : 
     508             :       TYPE(rho0_mpole_type), POINTER                     :: rho0_mpole
     509             :       CHARACTER(LEN=*), INTENT(IN)                       :: unit_str
     510             :       INTEGER, INTENT(in)                                :: output_unit
     511             : 
     512             :       INTEGER                                            :: ikind, l, nkind
     513             :       REAL(dp)                                           :: conv
     514             : 
     515           1 :       IF (ASSOCIATED(rho0_mpole)) THEN
     516           1 :          conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
     517             : 
     518             :          WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") &
     519           1 :             "*** Compensation density charges data set ***"
     520             :          WRITE (UNIT=output_unit, FMT="(T2,A,T35,f16.10)") &
     521           1 :             "- Rho0 exponent :", rho0_mpole%zet0_h
     522             :          WRITE (UNIT=output_unit, FMT="(T2,A,T35,I5)") &
     523           1 :             "- Global max l :", rho0_mpole%lmax_0
     524             : 
     525             :          WRITE (UNIT=output_unit, FMT="(T2,A)") &
     526           1 :             "- Normalization constants for g0"
     527           2 :          DO l = 0, rho0_mpole%lmax_0
     528             :             WRITE (UNIT=output_unit, FMT="(T20,A,T31,I2,T38,A,f15.5)") &
     529           2 :                "ang. mom.= ", l, " hard= ", rho0_mpole%norm_g0l_h(l)
     530             :          END DO
     531             : 
     532           1 :          nkind = SIZE(rho0_mpole%lmax0_kind, 1)
     533           2 :          DO ikind = 1, nkind
     534             :             WRITE (UNIT=output_unit, FMT="(/,T2,A,T55,I2)") &
     535             :                "- rho0 max L and radii in "//TRIM(unit_str)// &
     536           1 :                " for the atom kind :", ikind
     537             : 
     538             :             WRITE (UNIT=output_unit, FMT="(T2,T20,A,T55,I5)") &
     539           1 :                "=> l max  :", rho0_mpole%lmax0_kind(ikind)
     540             : 
     541             :             WRITE (UNIT=output_unit, FMT="(T2,T20,A,T55,f20.10)") &
     542           1 :                "=> max radius of g0: ", &
     543           3 :                rho0_mpole%mp_gau(ikind)%rpgf0_h*conv
     544             :          END DO ! ikind
     545             : 
     546             :       ELSE
     547             :          WRITE (UNIT=output_unit, FMT="(/,T5,A,/)") &
     548           0 :             ' WARNING: I cannot print rho0, it is not associated'
     549             :       END IF
     550             : 
     551           1 :    END SUBROUTINE write_rho0_info
     552           0 : END MODULE qs_rho0_types