Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : MODULE qs_kernel_methods
9 : USE cp_control_types, ONLY: dft_control_type
10 : USE cp_dbcsr_api, ONLY: dbcsr_distribution_type,&
11 : dbcsr_init_p,&
12 : dbcsr_p_type
13 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set
14 : USE input_section_types, ONLY: section_get_ival,&
15 : section_get_lval,&
16 : section_get_rval,&
17 : section_vals_get_subs_vals,&
18 : section_vals_type
19 : USE kinds, ONLY: dp
20 : USE pw_env_types, ONLY: pw_env_get,&
21 : pw_env_type
22 : USE pw_methods, ONLY: pw_axpy,&
23 : pw_zero
24 : USE pw_pool_types, ONLY: pw_pool_type
25 : USE pw_types, ONLY: pw_c1d_gs_type,&
26 : pw_r3d_rs_type
27 : USE qs_environment_types, ONLY: get_qs_env,&
28 : qs_environment_type
29 : USE qs_kernel_types, ONLY: full_kernel_env_type
30 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
31 : USE qs_rho_methods, ONLY: qs_rho_rebuild
32 : USE qs_rho_types, ONLY: qs_rho_create,&
33 : qs_rho_get,&
34 : qs_rho_set,&
35 : qs_rho_type
36 : USE qs_tddfpt2_subgroups, ONLY: tddfpt_dbcsr_create_by_dist,&
37 : tddfpt_subgroup_env_type
38 : USE xc, ONLY: xc_prep_2nd_deriv
39 : USE xc_derivatives, ONLY: xc_functionals_get_needs
40 : USE xc_fxc_kernel, ONLY: calc_fxc_kernel
41 : USE xc_rho_set_types, ONLY: xc_rho_set_create
42 : #include "./base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 :
46 : PRIVATE
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_kernel_methods'
49 :
50 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
51 :
52 : PUBLIC :: create_kernel_env, create_fxc_kernel
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief Create kernel environment.
58 : !> \param kernel_env kernel environment (allocated and initialised on exit)
59 : !> \param xc_section input section which defines an exchange-correlation functional
60 : !> \param is_rks_triplets indicates that the triplet excited states calculation using
61 : !> spin-unpolarised molecular orbitals has been requested
62 : !> \param rho_struct_sub ground state charge density, if not associated on input, it will be associated on output
63 : !> \param lsd_singlets ...
64 : !> \param do_excitations ...
65 : !> \param sub_env parallel group environment
66 : !> \param qs_env ...
67 : !> \par History
68 : !> * 02.2017 created [Sergey Chulkov]
69 : !> * 06.2018 the charge density needs to be provided via a dummy argument [Sergey Chulkov]
70 : ! **************************************************************************************************
71 988 : SUBROUTINE create_kernel_env(kernel_env, xc_section, is_rks_triplets, rho_struct_sub, &
72 : lsd_singlets, do_excitations, sub_env, qs_env)
73 : TYPE(full_kernel_env_type), INTENT(inout) :: kernel_env
74 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
75 : LOGICAL, INTENT(in) :: is_rks_triplets
76 : TYPE(qs_rho_type), POINTER :: rho_struct_sub
77 : LOGICAL, INTENT(in), OPTIONAL :: lsd_singlets, do_excitations
78 : TYPE(tddfpt_subgroup_env_type), INTENT(in), &
79 : OPTIONAL :: sub_env
80 : TYPE(qs_environment_type), INTENT(in), OPTIONAL, &
81 : POINTER :: qs_env
82 :
83 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_kernel_env'
84 :
85 : INTEGER :: handle, ispin, nspins
86 : LOGICAL :: lsd, my_excitations, my_singlets
87 : TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
88 988 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, rho_ia_ao
89 : TYPE(dft_control_type), POINTER :: dft_control
90 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
91 988 : POINTER :: sab_orb
92 : TYPE(pw_env_type), POINTER :: pw_env
93 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
94 988 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_ij_r, rho_ij_r2, tau_ij_r, tau_ij_r2
95 : TYPE(pw_r3d_rs_type), POINTER :: weights
96 : TYPE(section_vals_type), POINTER :: xc_fun_section
97 :
98 988 : CALL timeset(routineN, handle)
99 :
100 988 : IF (PRESENT(sub_env)) THEN
101 988 : pw_env => sub_env%pw_env
102 988 : dbcsr_dist => sub_env%dbcsr_dist
103 988 : sab_orb => sub_env%sab_orb
104 :
105 988 : nspins = SIZE(sub_env%mos_occ)
106 :
107 : NULLIFY (weights)
108 988 : weights => sub_env%xcint_weights
109 988 : IF (sub_env%is_split .AND. ASSOCIATED(weights)) THEN
110 0 : CPABORT("Split communicators and integration weights")
111 : END IF
112 :
113 : ELSE
114 0 : CPASSERT(PRESENT(qs_env))
115 :
116 : CALL get_qs_env(qs_env=qs_env, &
117 : pw_env=pw_env, &
118 : dbcsr_dist=dbcsr_dist, &
119 : sab_orb=sab_orb, &
120 0 : dft_control=dft_control)
121 :
122 0 : NULLIFY (weights)
123 0 : CALL get_qs_env(qs_env=qs_env, xcint_weights=weights)
124 :
125 0 : nspins = dft_control%nspins
126 :
127 0 : IF (.NOT. ASSOCIATED(rho_struct_sub)) THEN
128 : ! Build rho_set
129 0 : NULLIFY (rho_ia_ao)
130 0 : CALL dbcsr_allocate_matrix_set(rho_ia_ao, nspins)
131 0 : DO ispin = 1, nspins
132 0 : CALL dbcsr_init_p(rho_ia_ao(ispin)%matrix)
133 : CALL tddfpt_dbcsr_create_by_dist(rho_ia_ao(ispin)%matrix, template=matrix_s(1)%matrix, &
134 0 : dbcsr_dist=dbcsr_dist, sab=sab_orb)
135 : END DO
136 :
137 0 : ALLOCATE (rho_struct_sub)
138 0 : CALL qs_rho_create(rho_struct_sub)
139 0 : CALL qs_rho_set(rho_struct_sub, rho_ao=rho_ia_ao)
140 : CALL qs_rho_rebuild(rho_struct_sub, qs_env, rebuild_ao=.FALSE., &
141 0 : rebuild_grids=.TRUE., pw_env_external=pw_env)
142 : END IF
143 : END IF
144 :
145 988 : lsd = (nspins > 1) .OR. is_rks_triplets
146 :
147 988 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
148 :
149 988 : CALL qs_rho_get(rho_struct_sub, rho_r=rho_ij_r, tau_r=tau_ij_r)
150 :
151 988 : NULLIFY (kernel_env%xc_rho_set, kernel_env%xc_rho1_set)
152 :
153 22724 : ALLOCATE (kernel_env%xc_rho_set)
154 988 : IF (is_rks_triplets) THEN
155 : ! we are about to compute triplet states using spin-restricted reference MOs;
156 : ! we still need the beta-spin density component in order to compute the TDDFT kernel
157 330 : ALLOCATE (rho_ij_r2(2))
158 110 : rho_ij_r2(1) = rho_ij_r(1)
159 110 : rho_ij_r2(2) = rho_ij_r(1)
160 :
161 110 : IF (ASSOCIATED(tau_ij_r)) THEN
162 0 : ALLOCATE (tau_ij_r2(2))
163 0 : tau_ij_r2(1) = tau_ij_r(1)
164 0 : tau_ij_r2(2) = tau_ij_r(1)
165 : END IF
166 :
167 : CALL xc_prep_2nd_deriv(kernel_env%xc_deriv_set, kernel_env%xc_rho_set, rho_ij_r2, &
168 110 : auxbas_pw_pool, weights, xc_section=xc_section, tau_r=tau_ij_r2)
169 :
170 110 : IF (ASSOCIATED(tau_ij_r)) DEALLOCATE (tau_ij_r2)
171 :
172 110 : DEALLOCATE (rho_ij_r2)
173 : ELSE
174 : CALL xc_prep_2nd_deriv(kernel_env%xc_deriv_set, kernel_env%xc_rho_set, rho_ij_r, &
175 878 : auxbas_pw_pool, weights, xc_section=xc_section, tau_r=tau_ij_r)
176 : END IF
177 :
178 : ! ++ allocate structure for response density
179 988 : kernel_env%xc_section => xc_section
180 988 : kernel_env%deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV")
181 988 : kernel_env%rho_smooth_id = section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO")
182 :
183 988 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
184 : kernel_env%xc_rho1_cflags = xc_functionals_get_needs(functionals=xc_fun_section, lsd=lsd, &
185 988 : calc_potential=.TRUE.)
186 :
187 988 : IF (.NOT. ASSOCIATED(kernel_env%xc_rho1_set)) THEN
188 : NULLIFY (kernel_env%xc_rho1_set)
189 22724 : ALLOCATE (kernel_env%xc_rho1_set)
190 : END IF
191 : CALL xc_rho_set_create(kernel_env%xc_rho1_set, auxbas_pw_pool%pw_grid%bounds_local, &
192 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
193 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
194 988 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
195 :
196 988 : my_excitations = .TRUE.
197 988 : IF (PRESENT(do_excitations)) my_excitations = do_excitations
198 :
199 988 : my_singlets = .FALSE.
200 988 : IF (PRESENT(lsd_singlets)) my_singlets = lsd_singlets
201 :
202 988 : kernel_env%alpha = 1.0_dp
203 988 : kernel_env%beta = 0.0_dp
204 :
205 : ! kernel_env%beta is taken into account in spin-restricted case only
206 988 : IF (nspins == 1 .AND. my_excitations) THEN
207 864 : IF (is_rks_triplets) THEN
208 : ! K_{triplets} = K_{alpha,alpha} - K_{alpha,beta}
209 110 : kernel_env%beta = -1.0_dp
210 : ELSE
211 : ! alpha beta
212 : ! K_{singlets} = K_{alpha,alpha} + K_{alpha,beta} = 2 * K_{alpha,alpha} + 0 * K_{alpha,beta},
213 : ! due to the following relation : K_{alpha,alpha,singlets} == K_{alpha,beta,singlets}
214 754 : kernel_env%alpha = 2.0_dp
215 :
216 754 : IF (my_singlets) THEN
217 0 : kernel_env%beta = 1.0_dp
218 : END IF
219 : END IF
220 : END IF
221 :
222 : ! finite differences
223 988 : kernel_env%deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
224 988 : kernel_env%deriv3_analytic = section_get_lval(xc_section, "3RD_DERIV_ANALYTICAL")
225 :
226 988 : CALL timestop(handle)
227 :
228 988 : END SUBROUTINE create_kernel_env
229 :
230 : ! **************************************************************************************************
231 : !> \brief Create the xc kernel potential for the approximate Fxc kernel model
232 : !> \param rho_struct ...
233 : !> \param fxc_rspace ...
234 : !> \param xc_section ...
235 : !> \param is_rks_triplets ...
236 : !> \param sub_env ...
237 : !> \param qs_env ...
238 : ! **************************************************************************************************
239 24 : SUBROUTINE create_fxc_kernel(rho_struct, fxc_rspace, xc_section, is_rks_triplets, sub_env, qs_env)
240 : TYPE(qs_rho_type), POINTER :: rho_struct
241 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rspace
242 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
243 : LOGICAL, INTENT(IN) :: is_rks_triplets
244 : TYPE(tddfpt_subgroup_env_type), INTENT(IN) :: sub_env
245 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
246 :
247 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_fxc_kernel'
248 :
249 : INTEGER :: handle, ispin, nspins
250 : LOGICAL :: rho_g_valid, tau_r_valid
251 : REAL(KIND=dp) :: factor
252 12 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
253 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
254 : TYPE(pw_env_type), POINTER :: pw_env
255 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
256 12 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
257 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc
258 : TYPE(section_vals_type), POINTER :: xc_kernel
259 :
260 12 : CALL timeset(routineN, handle)
261 :
262 12 : pw_env => sub_env%pw_env
263 12 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
264 :
265 12 : NULLIFY (rho_r, rho_g, tau_r)
266 : CALL qs_rho_get(rho_struct, &
267 : tau_r_valid=tau_r_valid, &
268 : rho_g_valid=rho_g_valid, &
269 : rho_r=rho_r, &
270 : rho_g=rho_g, &
271 12 : tau_r=tau_r)
272 :
273 12 : IF (.NOT. tau_r_valid) NULLIFY (tau_r)
274 12 : IF (.NOT. rho_g_valid) NULLIFY (rho_g)
275 :
276 12 : nspins = SIZE(rho_r)
277 :
278 12 : NULLIFY (rho_nlcc, rho_nlcc_g)
279 : CALL get_qs_env(qs_env, &
280 : rho_nlcc=rho_nlcc, &
281 12 : rho_nlcc_g=rho_nlcc_g)
282 : ! add the nlcc densities
283 12 : IF (ASSOCIATED(rho_nlcc)) THEN
284 0 : factor = 1.0_dp
285 0 : DO ispin = 1, nspins
286 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
287 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
288 : END DO
289 : END IF
290 :
291 48 : DO ispin = 1, SIZE(fxc_rspace)
292 48 : CALL pw_zero(fxc_rspace(ispin))
293 : END DO
294 :
295 12 : xc_kernel => section_vals_get_subs_vals(xc_section, "XC_KERNEL")
296 : CALL calc_fxc_kernel(fxc_rspace, rho_r, rho_g, tau_r, &
297 12 : xc_kernel, is_rks_triplets, auxbas_pw_pool)
298 :
299 : ! remove the nlcc densities (keep stuff in original state)
300 12 : IF (ASSOCIATED(rho_nlcc)) THEN
301 0 : factor = -1.0_dp
302 0 : DO ispin = 1, nspins
303 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
304 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
305 : END DO
306 : END IF
307 :
308 12 : CALL timestop(handle)
309 :
310 12 : END SUBROUTINE create_fxc_kernel
311 :
312 : END MODULE qs_kernel_methods
|
