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 628 : 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 628 : 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 628 : CALL timeset(routineN, handle)
87 :
88 628 : nspins = SIZE(sub_env%mos_occ)
89 1376 : DO ispin = 1, nspins
90 1376 : CALL cp_fm_get_info(sub_env%mos_occ(ispin), nrow_global=nao, ncol_global=nmo_occ(ispin))
91 : END DO
92 :
93 628 : CALL qs_rho_get(rho_orb_struct, rho_ao=rho_ij_ao)
94 1376 : 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 748 : 0.0_dp, wfm_rho_orb)
98 :
99 1376 : 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 628 : IF (nspins == 1 .AND. (.NOT. is_rks_triplets)) &
104 428 : CALL dbcsr_scale(rho_ij_ao(1)%matrix, 2.0_dp)
105 :
106 628 : CALL get_qs_env(qs_env, dft_control=dft_control)
107 :
108 628 : 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 188 : para_env_external=sub_env%para_env)
114 188 : CALL prepare_gapw_den(qs_env, local_rho_set=sub_env%local_rho_set, pw_env_sub=sub_env%pw_env)
115 440 : 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 404 : para_env_external=sub_env%para_env)
132 : END IF
133 :
134 628 : CALL timestop(handle)
135 :
136 628 : 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 1948 : 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 974 : REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_aux_fit_r
173 : TYPE(admm_type), POINTER :: admm_env
174 974 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao_aux_fit, rho_ao_orb
175 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
176 974 : POINTER :: sab_aux_fit
177 974 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_aux_fit_g
178 974 : 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 974 : CALL timeset(routineN, handle)
183 :
184 974 : CPASSERT(ASSOCIATED(sub_env%admm_A))
185 :
186 974 : CALL get_qs_env(qs_env, ks_env=ks_env, admm_env=admm_env)
187 974 : 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 974 : rho_r=rho_aux_fit_r, tot_rho_r=tot_rho_aux_fit_r)
190 :
191 974 : nspins = SIZE(rho_ao_orb)
192 :
193 974 : IF (admm_env%do_gapw) THEN
194 176 : basis_type = "AUX_FIT_SOFT"
195 176 : task_list => sub_env%task_list_aux_fit_soft
196 : ELSE
197 798 : basis_type = "AUX_FIT"
198 798 : task_list => sub_env%task_list_aux_fit
199 : END IF
200 :
201 974 : CALL cp_fm_get_info(sub_env%admm_A, nrow_global=nao_aux, ncol_global=nao)
202 1952 : DO ispin = 1, nspins
203 : ! TO DO: consider sub_env%admm_A to be a DBCSR matrix
204 978 : 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 978 : 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 978 : 0.0_dp, wfm_rho_aux_fit)
209 978 : 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 1952 : pw_env_external=sub_env%pw_env, task_list_external=task_list)
216 : END DO
217 974 : IF (admm_env%do_gapw) THEN
218 176 : 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 176 : 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 176 : pw_env_sub=sub_env%pw_env)
226 : END IF
227 :
228 974 : CALL timestop(handle)
229 :
230 974 : END SUBROUTINE tddfpt_construct_aux_fit_density
231 :
232 : END MODULE qs_tddfpt2_densities
|
