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             : MODULE qs_tddfpt2_densities
       9             :    USE admm_types,                      ONLY: admm_type,&
      10             :                                               get_admm_env
      11             :    USE cp_control_types,                ONLY: dft_control_type
      12             :    USE cp_dbcsr_api,                    ONLY: dbcsr_p_type,&
      13             :                                               dbcsr_scale
      14             :    USE cp_dbcsr_operations,             ONLY: copy_dbcsr_to_fm,&
      15             :                                               copy_fm_to_dbcsr
      16             :    USE cp_fm_types,                     ONLY: cp_fm_get_info,&
      17             :                                               cp_fm_type
      18             :    USE kinds,                           ONLY: default_string_length,&
      19             :                                               dp
      20             :    USE parallel_gemm_api,               ONLY: parallel_gemm
      21             :    USE pw_env_types,                    ONLY: pw_env_get
      22             :    USE pw_pool_types,                   ONLY: pw_pool_type
      23             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      24             :                                               pw_r3d_rs_type
      25             :    USE qs_collocate_density,            ONLY: calculate_rho_elec
      26             :    USE qs_environment_types,            ONLY: get_qs_env,&
      27             :                                               qs_environment_type
      28             :    USE qs_gapw_densities,               ONLY: prepare_gapw_den
      29             :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      30             :    USE qs_local_rho_types,              ONLY: local_rho_type
      31             :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      32             :    USE qs_rho_atom_methods,             ONLY: calculate_rho_atom_coeff
      33             :    USE qs_rho_methods,                  ONLY: qs_rho_copy,&
      34             :                                               qs_rho_update_rho
      35             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      36             :                                               qs_rho_type
      37             :    USE qs_tddfpt2_subgroups,            ONLY: tddfpt_subgroup_env_type
      38             :    USE task_list_types,                 ONLY: task_list_type
      39             : #include "./base/base_uses.f90"
      40             : 
      41             :    IMPLICIT NONE
      42             : 
      43             :    PRIVATE
      44             : 
      45             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_tddfpt2_densities'
      46             : 
      47             :    INTEGER, PARAMETER, PRIVATE          :: maxspins = 2
      48             : 
      49             :    PUBLIC :: tddfpt_construct_ground_state_orb_density, tddfpt_construct_aux_fit_density
      50             : 
      51             : ! **************************************************************************************************
      52             : 
      53             : CONTAINS
      54             : 
      55             : ! **************************************************************************************************
      56             : !> \brief Compute the ground-state charge density expressed in primary basis set.
      57             : !> \param rho_orb_struct      ground-state density in primary basis set
      58             : !> \param rho_xc_struct       ground-state density in primary basis set for GAPW_XC
      59             : !> \param is_rks_triplets     indicates that the triplet excited states calculation using
      60             : !>                            spin-unpolarised molecular orbitals has been requested
      61             : !> \param qs_env              Quickstep environment
      62             : !> \param sub_env             parallel (sub)group environment
      63             : !> \param wfm_rho_orb         work dense matrix with shape [nao x nao] distributed among
      64             : !>                            processors of the given parallel group (modified on exit)
      65             : !> \par History
      66             : !>    * 06.2018 created by splitting the subroutine tddfpt_apply_admm_correction() in two
      67             : !>              subroutines tddfpt_construct_ground_state_orb_density() and
      68             : !>              tddfpt_construct_aux_fit_density [Sergey Chulkov]
      69             : ! **************************************************************************************************
      70         578 :    SUBROUTINE tddfpt_construct_ground_state_orb_density(rho_orb_struct, rho_xc_struct, is_rks_triplets, &
      71             :                                                         qs_env, sub_env, wfm_rho_orb)
      72             :       TYPE(qs_rho_type), POINTER                         :: rho_orb_struct, rho_xc_struct
      73             :       LOGICAL, INTENT(in)                                :: is_rks_triplets
      74             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      75             :       TYPE(tddfpt_subgroup_env_type), INTENT(in)         :: sub_env
      76             :       TYPE(cp_fm_type), INTENT(IN)                       :: wfm_rho_orb
      77             : 
      78             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_construct_ground_state_orb_density'
      79             : 
      80             :       INTEGER                                            :: handle, ispin, nao, nspins
      81             :       INTEGER, DIMENSION(maxspins)                       :: nmo_occ
      82         578 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ij_ao
      83             :       TYPE(dft_control_type), POINTER                    :: dft_control
      84             :       TYPE(pw_pool_type), POINTER                        :: auxbas_pw_pool
      85             : 
      86         578 :       CALL timeset(routineN, handle)
      87             : 
      88         578 :       nspins = SIZE(sub_env%mos_occ)
      89        1248 :       DO ispin = 1, nspins
      90        1248 :          CALL cp_fm_get_info(sub_env%mos_occ(ispin), nrow_global=nao, ncol_global=nmo_occ(ispin))
      91             :       END DO
      92             : 
      93         578 :       CALL qs_rho_get(rho_orb_struct, rho_ao=rho_ij_ao)
      94        1248 :       DO ispin = 1, nspins
      95             :          CALL parallel_gemm('N', 'T', nao, nao, nmo_occ(ispin), 1.0_dp, &
      96             :                             sub_env%mos_occ(ispin), sub_env%mos_occ(ispin), &
      97         670 :                             0.0_dp, wfm_rho_orb)
      98             : 
      99        1248 :          CALL copy_fm_to_dbcsr(wfm_rho_orb, rho_ij_ao(ispin)%matrix, keep_sparsity=.TRUE.)
     100             :       END DO
     101             : 
     102             :       ! take into account that all MOs are doubly occupied in spin-restricted case
     103         578 :       IF (nspins == 1 .AND. (.NOT. is_rks_triplets)) &
     104         406 :          CALL dbcsr_scale(rho_ij_ao(1)%matrix, 2.0_dp)
     105             : 
     106         578 :       CALL get_qs_env(qs_env, dft_control=dft_control)
     107             : 
     108         578 :       IF (dft_control%qs_control%gapw) THEN
     109             :          CALL qs_rho_update_rho(rho_orb_struct, qs_env, &
     110             :                                 local_rho_set=sub_env%local_rho_set, &
     111             :                                 pw_env_external=sub_env%pw_env, &
     112             :                                 task_list_external=sub_env%task_list_orb_soft, &
     113         164 :                                 para_env_external=sub_env%para_env)
     114         164 :          CALL prepare_gapw_den(qs_env, local_rho_set=sub_env%local_rho_set, pw_env_sub=sub_env%pw_env)
     115         414 :       ELSEIF (dft_control%qs_control%gapw_xc) THEN
     116             :          CALL qs_rho_update_rho(rho_orb_struct, qs_env, &
     117             :                                 rho_xc_external=rho_xc_struct, &
     118             :                                 local_rho_set=sub_env%local_rho_set, &
     119             :                                 pw_env_external=sub_env%pw_env, &
     120             :                                 task_list_external=sub_env%task_list_orb, &
     121             :                                 task_list_external_soft=sub_env%task_list_orb_soft, &
     122          36 :                                 para_env_external=sub_env%para_env)
     123          36 :          CALL pw_env_get(sub_env%pw_env, auxbas_pw_pool=auxbas_pw_pool)
     124          36 :          CALL qs_rho_copy(rho_xc_struct, rho_orb_struct, auxbas_pw_pool, nspins)
     125             :          CALL prepare_gapw_den(qs_env, local_rho_set=sub_env%local_rho_set, do_rho0=.FALSE., &
     126          36 :                                pw_env_sub=sub_env%pw_env)
     127             :       ELSE
     128             :          CALL qs_rho_update_rho(rho_orb_struct, qs_env, &
     129             :                                 pw_env_external=sub_env%pw_env, &
     130             :                                 task_list_external=sub_env%task_list_orb, &
     131         378 :                                 para_env_external=sub_env%para_env)
     132             :       END IF
     133             : 
     134         578 :       CALL timestop(handle)
     135             : 
     136         578 :    END SUBROUTINE tddfpt_construct_ground_state_orb_density
     137             : 
     138             : ! **************************************************************************************************
     139             : !> \brief Project a charge density expressed in primary basis set into the auxiliary basis set.
     140             : !> \param rho_orb_struct      response density in primary basis set
     141             : !> \param rho_aux_fit_struct  response density in auxiliary basis set (modified on exit)
     142             : !> \param local_rho_set       GAPW density of auxiliary basis set density
     143             : !> \param qs_env              Quickstep environment
     144             : !> \param sub_env             parallel (sub)group environment
     145             : !> \param wfm_rho_orb         work dense matrix with shape [nao x nao] distributed among
     146             : !>                            processors of the given parallel group (modified on exit)
     147             : !> \param wfm_rho_aux_fit     work dense matrix with shape [nao_aux x nao_aux] distributed among
     148             : !>                            processors of the given parallel group (modified on exit)
     149             : !> \param wfm_aux_orb         work dense matrix with shape [nao_aux x nao] distributed among
     150             : !>                            processors of the given parallel group (modified on exit)
     151             : !> \par History
     152             : !>    * 03.2017 the subroutine tddfpt_apply_admm_correction() was originally created by splitting
     153             : !>              the subroutine tddfpt_apply_hfx() in two parts [Sergey Chulkov]
     154             : !>    * 06.2018 created by splitting the subroutine tddfpt_apply_admm_correction() in two subroutines
     155             : !>              tddfpt_construct_ground_state_orb_density() and tddfpt_construct_aux_fit_density()
     156             : !>              in order to avoid code duplication [Sergey Chulkov]
     157             : ! **************************************************************************************************
     158        1816 :    SUBROUTINE tddfpt_construct_aux_fit_density(rho_orb_struct, rho_aux_fit_struct, local_rho_set, &
     159             :                                                qs_env, sub_env, &
     160             :                                                wfm_rho_orb, wfm_rho_aux_fit, wfm_aux_orb)
     161             :       TYPE(qs_rho_type), POINTER                         :: rho_orb_struct, rho_aux_fit_struct
     162             :       TYPE(local_rho_type), POINTER                      :: local_rho_set
     163             :       TYPE(qs_environment_type), POINTER                 :: qs_env
     164             :       TYPE(tddfpt_subgroup_env_type), INTENT(in)         :: sub_env
     165             :       TYPE(cp_fm_type), INTENT(INOUT)                    :: wfm_rho_orb
     166             :       TYPE(cp_fm_type), INTENT(IN)                       :: wfm_rho_aux_fit, wfm_aux_orb
     167             : 
     168             :       CHARACTER(LEN=*), PARAMETER :: routineN = 'tddfpt_construct_aux_fit_density'
     169             : 
     170             :       CHARACTER(LEN=default_string_length)               :: basis_type
     171             :       INTEGER                                            :: handle, ispin, nao, nao_aux, nspins
     172         908 :       REAL(kind=dp), DIMENSION(:), POINTER               :: tot_rho_aux_fit_r
     173             :       TYPE(admm_type), POINTER                           :: admm_env
     174         908 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: rho_ao_aux_fit, rho_ao_orb
     175             :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     176         908 :          POINTER                                         :: sab_aux_fit
     177         908 :       TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER        :: rho_aux_fit_g
     178         908 :       TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER        :: rho_aux_fit_r
     179             :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     180             :       TYPE(task_list_type), POINTER                      :: task_list
     181             : 
     182         908 :       CALL timeset(routineN, handle)
     183             : 
     184         908 :       CPASSERT(ASSOCIATED(sub_env%admm_A))
     185             : 
     186         908 :       CALL get_qs_env(qs_env, ks_env=ks_env, admm_env=admm_env)
     187         908 :       CALL qs_rho_get(rho_orb_struct, rho_ao=rho_ao_orb)
     188             :       CALL qs_rho_get(rho_aux_fit_struct, rho_ao=rho_ao_aux_fit, rho_g=rho_aux_fit_g, &
     189         908 :                       rho_r=rho_aux_fit_r, tot_rho_r=tot_rho_aux_fit_r)
     190             : 
     191         908 :       nspins = SIZE(rho_ao_orb)
     192             : 
     193         908 :       IF (admm_env%do_gapw) THEN
     194         136 :          basis_type = "AUX_FIT_SOFT"
     195         136 :          task_list => sub_env%task_list_aux_fit_soft
     196             :       ELSE
     197         772 :          basis_type = "AUX_FIT"
     198         772 :          task_list => sub_env%task_list_aux_fit
     199             :       END IF
     200             : 
     201         908 :       CALL cp_fm_get_info(sub_env%admm_A, nrow_global=nao_aux, ncol_global=nao)
     202        1820 :       DO ispin = 1, nspins
     203             :          ! TO DO: consider sub_env%admm_A to be a DBCSR matrix
     204         912 :          CALL copy_dbcsr_to_fm(rho_ao_orb(ispin)%matrix, wfm_rho_orb)
     205             :          CALL parallel_gemm('N', 'N', nao_aux, nao, nao, 1.0_dp, sub_env%admm_A, &
     206         912 :                             wfm_rho_orb, 0.0_dp, wfm_aux_orb)
     207             :          CALL parallel_gemm('N', 'T', nao_aux, nao_aux, nao, 1.0_dp, sub_env%admm_A, wfm_aux_orb, &
     208         912 :                             0.0_dp, wfm_rho_aux_fit)
     209         912 :          CALL copy_fm_to_dbcsr(wfm_rho_aux_fit, rho_ao_aux_fit(ispin)%matrix, keep_sparsity=.TRUE.)
     210             : 
     211             :          CALL calculate_rho_elec(matrix_p=rho_ao_aux_fit(ispin)%matrix, &
     212             :                                  rho=rho_aux_fit_r(ispin), rho_gspace=rho_aux_fit_g(ispin), &
     213             :                                  total_rho=tot_rho_aux_fit_r(ispin), ks_env=ks_env, &
     214             :                                  soft_valid=.FALSE., basis_type=basis_type, &
     215        1820 :                                  pw_env_external=sub_env%pw_env, task_list_external=task_list)
     216             :       END DO
     217         908 :       IF (admm_env%do_gapw) THEN
     218         136 :          CALL get_admm_env(qs_env%admm_env, sab_aux_fit=sab_aux_fit)
     219             :          CALL calculate_rho_atom_coeff(qs_env, rho_ao_aux_fit, &
     220             :                                        rho_atom_set=local_rho_set%rho_atom_set, &
     221             :                                        qs_kind_set=admm_env%admm_gapw_env%admm_kind_set, &
     222         136 :                                        oce=admm_env%admm_gapw_env%oce, sab=sab_aux_fit, para_env=sub_env%para_env)
     223             :          CALL prepare_gapw_den(qs_env, local_rho_set=local_rho_set, &
     224             :                                do_rho0=.FALSE., kind_set_external=admm_env%admm_gapw_env%admm_kind_set, &
     225         136 :                                pw_env_sub=sub_env%pw_env)
     226             :       END IF
     227             : 
     228         908 :       CALL timestop(handle)
     229             : 
     230         908 :    END SUBROUTINE tddfpt_construct_aux_fit_density
     231             : 
     232             : END MODULE qs_tddfpt2_densities