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_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 756 : 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 756 : 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 756 : POINTER :: sab_orb
92 : TYPE(pw_env_type), POINTER :: pw_env
93 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
94 756 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_ij_r, rho_ij_r2, tau_ij_r, tau_ij_r2
95 : TYPE(section_vals_type), POINTER :: xc_fun_section
96 :
97 756 : CALL timeset(routineN, handle)
98 :
99 756 : IF (PRESENT(sub_env)) THEN
100 756 : pw_env => sub_env%pw_env
101 756 : dbcsr_dist => sub_env%dbcsr_dist
102 756 : sab_orb => sub_env%sab_orb
103 :
104 756 : nspins = SIZE(sub_env%mos_occ)
105 : ELSE
106 0 : CPASSERT(PRESENT(qs_env))
107 :
108 : CALL get_qs_env(qs_env=qs_env, &
109 : pw_env=pw_env, &
110 : dbcsr_dist=dbcsr_dist, &
111 : sab_orb=sab_orb, &
112 0 : dft_control=dft_control)
113 :
114 0 : nspins = dft_control%nspins
115 :
116 0 : IF (.NOT. ASSOCIATED(rho_struct_sub)) THEN
117 : ! Build rho_set
118 0 : NULLIFY (rho_ia_ao)
119 0 : CALL dbcsr_allocate_matrix_set(rho_ia_ao, nspins)
120 0 : DO ispin = 1, nspins
121 0 : CALL dbcsr_init_p(rho_ia_ao(ispin)%matrix)
122 : CALL tddfpt_dbcsr_create_by_dist(rho_ia_ao(ispin)%matrix, template=matrix_s(1)%matrix, &
123 0 : dbcsr_dist=dbcsr_dist, sab=sab_orb)
124 : END DO
125 :
126 0 : ALLOCATE (rho_struct_sub)
127 0 : CALL qs_rho_create(rho_struct_sub)
128 0 : CALL qs_rho_set(rho_struct_sub, rho_ao=rho_ia_ao)
129 : CALL qs_rho_rebuild(rho_struct_sub, qs_env, rebuild_ao=.FALSE., &
130 0 : rebuild_grids=.TRUE., pw_env_external=pw_env)
131 : END IF
132 : END IF
133 :
134 756 : lsd = (nspins > 1) .OR. is_rks_triplets
135 :
136 756 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
137 :
138 756 : CALL qs_rho_get(rho_struct_sub, rho_r=rho_ij_r, tau_r=tau_ij_r)
139 :
140 756 : NULLIFY (kernel_env%xc_rho_set, kernel_env%xc_rho1_set)
141 :
142 17388 : ALLOCATE (kernel_env%xc_rho_set)
143 756 : IF (is_rks_triplets) THEN
144 : ! we are about to compute triplet states using spin-restricted reference MOs;
145 : ! we still need the beta-spin density component in order to compute the TDDFT kernel
146 240 : ALLOCATE (rho_ij_r2(2))
147 80 : rho_ij_r2(1) = rho_ij_r(1)
148 80 : rho_ij_r2(2) = rho_ij_r(1)
149 :
150 80 : IF (ASSOCIATED(tau_ij_r)) THEN
151 0 : ALLOCATE (tau_ij_r2(2))
152 0 : tau_ij_r2(1) = tau_ij_r(1)
153 0 : tau_ij_r2(2) = tau_ij_r(1)
154 : END IF
155 :
156 : CALL xc_prep_2nd_deriv(kernel_env%xc_deriv_set, kernel_env%xc_rho_set, rho_ij_r2, &
157 80 : auxbas_pw_pool, xc_section=xc_section, tau_r=tau_ij_r2)
158 :
159 80 : IF (ASSOCIATED(tau_ij_r)) DEALLOCATE (tau_ij_r2)
160 :
161 80 : DEALLOCATE (rho_ij_r2)
162 : ELSE
163 : CALL xc_prep_2nd_deriv(kernel_env%xc_deriv_set, kernel_env%xc_rho_set, rho_ij_r, &
164 676 : auxbas_pw_pool, xc_section=xc_section, tau_r=tau_ij_r)
165 : END IF
166 :
167 : ! ++ allocate structure for response density
168 756 : kernel_env%xc_section => xc_section
169 756 : kernel_env%deriv_method_id = section_get_ival(xc_section, "XC_GRID%XC_DERIV")
170 756 : kernel_env%rho_smooth_id = section_get_ival(xc_section, "XC_GRID%XC_SMOOTH_RHO")
171 :
172 756 : xc_fun_section => section_vals_get_subs_vals(xc_section, "XC_FUNCTIONAL")
173 : kernel_env%xc_rho1_cflags = xc_functionals_get_needs(functionals=xc_fun_section, lsd=lsd, &
174 756 : calc_potential=.TRUE.)
175 :
176 756 : IF (.NOT. ASSOCIATED(kernel_env%xc_rho1_set)) THEN
177 : NULLIFY (kernel_env%xc_rho1_set)
178 17388 : ALLOCATE (kernel_env%xc_rho1_set)
179 : END IF
180 : CALL xc_rho_set_create(kernel_env%xc_rho1_set, auxbas_pw_pool%pw_grid%bounds_local, &
181 : rho_cutoff=section_get_rval(xc_section, "DENSITY_CUTOFF"), &
182 : drho_cutoff=section_get_rval(xc_section, "GRADIENT_CUTOFF"), &
183 756 : tau_cutoff=section_get_rval(xc_section, "TAU_CUTOFF"))
184 :
185 756 : my_excitations = .TRUE.
186 756 : IF (PRESENT(do_excitations)) my_excitations = do_excitations
187 :
188 756 : my_singlets = .FALSE.
189 756 : IF (PRESENT(lsd_singlets)) my_singlets = lsd_singlets
190 :
191 756 : kernel_env%alpha = 1.0_dp
192 756 : kernel_env%beta = 0.0_dp
193 :
194 : ! kernel_env%beta is taken into account in spin-restricted case only
195 756 : IF (nspins == 1 .AND. my_excitations) THEN
196 632 : IF (is_rks_triplets) THEN
197 : ! K_{triplets} = K_{alpha,alpha} - K_{alpha,beta}
198 80 : kernel_env%beta = -1.0_dp
199 : ELSE
200 : ! alpha beta
201 : ! K_{singlets} = K_{alpha,alpha} + K_{alpha,beta} = 2 * K_{alpha,alpha} + 0 * K_{alpha,beta},
202 : ! due to the following relation : K_{alpha,alpha,singlets} == K_{alpha,beta,singlets}
203 552 : kernel_env%alpha = 2.0_dp
204 :
205 552 : IF (my_singlets) THEN
206 0 : kernel_env%beta = 1.0_dp
207 : END IF
208 : END IF
209 : END IF
210 :
211 : ! finite differences
212 756 : kernel_env%deriv2_analytic = section_get_lval(xc_section, "2ND_DERIV_ANALYTICAL")
213 756 : kernel_env%deriv3_analytic = section_get_lval(xc_section, "3RD_DERIV_ANALYTICAL")
214 :
215 756 : CALL timestop(handle)
216 :
217 756 : END SUBROUTINE create_kernel_env
218 :
219 : ! **************************************************************************************************
220 : !> \brief Create the xc kernel potential for the approximate Fxc kernel model
221 : !> \param rho_struct ...
222 : !> \param fxc_rspace ...
223 : !> \param xc_section ...
224 : !> \param is_rks_triplets ...
225 : !> \param sub_env ...
226 : !> \param qs_env ...
227 : ! **************************************************************************************************
228 24 : SUBROUTINE create_fxc_kernel(rho_struct, fxc_rspace, xc_section, is_rks_triplets, sub_env, qs_env)
229 : TYPE(qs_rho_type), POINTER :: rho_struct
230 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fxc_rspace
231 : TYPE(section_vals_type), INTENT(IN), POINTER :: xc_section
232 : LOGICAL, INTENT(IN) :: is_rks_triplets
233 : TYPE(tddfpt_subgroup_env_type), INTENT(IN) :: sub_env
234 : TYPE(qs_environment_type), INTENT(IN), POINTER :: qs_env
235 :
236 : CHARACTER(LEN=*), PARAMETER :: routineN = 'create_fxc_kernel'
237 :
238 : INTEGER :: handle, ispin, nspins
239 : LOGICAL :: rho_g_valid, tau_r_valid
240 : REAL(KIND=dp) :: factor
241 12 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
242 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
243 : TYPE(pw_env_type), POINTER :: pw_env
244 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
245 12 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, tau_r
246 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc
247 : TYPE(section_vals_type), POINTER :: xc_kernel
248 :
249 12 : CALL timeset(routineN, handle)
250 :
251 12 : pw_env => sub_env%pw_env
252 12 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
253 :
254 12 : NULLIFY (rho_r, rho_g, tau_r)
255 : CALL qs_rho_get(rho_struct, &
256 : tau_r_valid=tau_r_valid, &
257 : rho_g_valid=rho_g_valid, &
258 : rho_r=rho_r, &
259 : rho_g=rho_g, &
260 12 : tau_r=tau_r)
261 :
262 12 : IF (.NOT. tau_r_valid) NULLIFY (tau_r)
263 12 : IF (.NOT. rho_g_valid) NULLIFY (rho_g)
264 :
265 12 : nspins = SIZE(rho_r)
266 :
267 12 : NULLIFY (rho_nlcc, rho_nlcc_g)
268 : CALL get_qs_env(qs_env, &
269 : rho_nlcc=rho_nlcc, &
270 12 : rho_nlcc_g=rho_nlcc_g)
271 : ! add the nlcc densities
272 12 : IF (ASSOCIATED(rho_nlcc)) THEN
273 0 : factor = 1.0_dp
274 0 : DO ispin = 1, nspins
275 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
276 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
277 : END DO
278 : END IF
279 :
280 48 : DO ispin = 1, SIZE(fxc_rspace)
281 48 : CALL pw_zero(fxc_rspace(ispin))
282 : END DO
283 :
284 12 : xc_kernel => section_vals_get_subs_vals(xc_section, "XC_KERNEL")
285 : CALL calc_fxc_kernel(fxc_rspace, rho_r, rho_g, tau_r, &
286 12 : xc_kernel, is_rks_triplets, auxbas_pw_pool)
287 :
288 : ! remove the nlcc densities (keep stuff in original state)
289 12 : IF (ASSOCIATED(rho_nlcc)) THEN
290 0 : factor = -1.0_dp
291 0 : DO ispin = 1, nspins
292 0 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
293 0 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
294 : END DO
295 : END IF
296 :
297 12 : CALL timestop(handle)
298 :
299 12 : END SUBROUTINE create_fxc_kernel
300 :
301 : END MODULE qs_kernel_methods
|
