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 : ! **************************************************************************************************
9 : !> \brief Calculate the derivatives of the MO coefficients wrt nuclear coordinates
10 : !> \author Sandra Luber, Edward Ditler
11 : ! **************************************************************************************************
12 :
13 : MODULE qs_dcdr_ao
14 :
15 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
16 : gto_basis_set_type
17 : USE cp_control_types, ONLY: dft_control_type
18 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
19 : dbcsr_get_block_p,&
20 : dbcsr_p_type,&
21 : dbcsr_set,&
22 : dbcsr_type
23 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
24 : copy_fm_to_dbcsr
25 : USE cp_fm_types, ONLY: cp_fm_create,&
26 : cp_fm_release,&
27 : cp_fm_type
28 : USE cp_log_handling, ONLY: cp_get_default_logger,&
29 : cp_logger_type
30 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
31 : section_vals_type
32 : USE kinds, ONLY: default_string_length,&
33 : dp
34 : USE orbital_pointers, ONLY: ncoset
35 : USE parallel_gemm_api, ONLY: parallel_gemm
36 : USE pw_env_types, ONLY: pw_env_get,&
37 : pw_env_type
38 : USE pw_methods, ONLY: pw_axpy,&
39 : pw_copy,&
40 : pw_scale,&
41 : pw_transfer,&
42 : pw_zero
43 : USE pw_poisson_methods, ONLY: pw_poisson_solve
44 : USE pw_poisson_types, ONLY: pw_poisson_type
45 : USE pw_pool_types, ONLY: pw_pool_p_type,&
46 : pw_pool_type
47 : USE pw_types, ONLY: pw_c1d_gs_type,&
48 : pw_r3d_rs_type
49 : USE qs_collocate_density, ONLY: calculate_drho_core,&
50 : calculate_drho_elec_dR
51 : USE qs_core_matrices, ONLY: core_matrices
52 : USE qs_energy_types, ONLY: qs_energy_type
53 : USE qs_environment_types, ONLY: get_qs_env,&
54 : qs_environment_type
55 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
56 : get_memory_usage
57 : USE qs_integrate_potential, ONLY: integrate_v_dbasis,&
58 : integrate_v_rspace
59 : USE qs_kind_types, ONLY: qs_kind_type
60 : USE qs_ks_types, ONLY: get_ks_env,&
61 : qs_ks_env_type
62 : USE qs_linres_types, ONLY: dcdr_env_type
63 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
64 : get_neighbor_list_set_p,&
65 : neighbor_list_iterate,&
66 : neighbor_list_iterator_create,&
67 : neighbor_list_iterator_p_type,&
68 : neighbor_list_iterator_release,&
69 : neighbor_list_set_p_type
70 : USE qs_rho_methods, ONLY: qs_rho_rebuild,&
71 : qs_rho_update_rho
72 : USE qs_rho_types, ONLY: qs_rho_create,&
73 : qs_rho_get,&
74 : qs_rho_release,&
75 : qs_rho_type
76 : USE qs_vxc, ONLY: qs_vxc_create
77 : USE xc, ONLY: xc_calc_2nd_deriv,&
78 : xc_prep_2nd_deriv
79 : USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
80 : xc_dset_release
81 : USE xc_rho_set_types, ONLY: xc_rho_set_release,&
82 : xc_rho_set_type
83 :
84 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
85 : !$ USE OMP_LIB, ONLY: omp_lock_kind, &
86 : !$ omp_init_lock, omp_set_lock, &
87 : !$ omp_unset_lock, omp_destroy_lock
88 :
89 : #include "./base/base_uses.f90"
90 :
91 : IMPLICIT NONE
92 :
93 : PRIVATE
94 : PUBLIC :: core_dR, d_vhxc_dR, d_core_charge_density_dR, apply_op_constant_term
95 : PUBLIC :: vhxc_R_perturbed_basis_functions
96 : PUBLIC :: hr_mult_by_delta_1d
97 :
98 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr_ao'
99 : CHARACTER(len=*), PARAMETER, PRIVATE :: dcdr_meta_gga_error = &
100 : "Analytical DCDR is not implemented for functionals that depend on the kinetic energy density. "// &
101 : "Use PROPERTIES%LINRES%DCDR%APT_FD T to calculate APTs by finite differences."
102 :
103 : CONTAINS
104 :
105 : ! **************************************************************************************************
106 : !> \brief Build the perturbed density matrix correction depending on the overlap derivative
107 : !> \param qs_env ...
108 : !> \param dcdr_env ...
109 : !> \param overlap1 Overlap derivative in AO basis
110 : !> \author Edward Ditler
111 : ! **************************************************************************************************
112 252 : SUBROUTINE apply_op_constant_term(qs_env, dcdr_env, overlap1)
113 : TYPE(qs_environment_type), POINTER :: qs_env
114 : TYPE(dcdr_env_type) :: dcdr_env
115 : TYPE(dbcsr_p_type), OPTIONAL :: overlap1
116 :
117 : CHARACTER(len=*), PARAMETER :: routineN = 'apply_op_constant_term'
118 :
119 : INTEGER :: handle, ispin
120 : REAL(KIND=dp) :: energy_hartree
121 : TYPE(cp_fm_type) :: rho_ao_fm, rho_ao_s1, rho_ao_s1_rho_ao, &
122 : s1_ao
123 252 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho1_ao, rho_ao
124 : TYPE(pw_c1d_gs_type) :: rho1_tot_gspace, v_hartree_gspace
125 252 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g, rho1_g_pw
126 : TYPE(pw_env_type), POINTER :: pw_env
127 : TYPE(pw_poisson_type), POINTER :: poisson_env
128 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
129 : TYPE(pw_r3d_rs_type) :: v_hartree_rspace
130 504 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho_r, tau1_r, v_rspace_new, &
131 252 : v_xc, v_xc_tau
132 : TYPE(pw_r3d_rs_type), POINTER :: weights
133 : TYPE(qs_rho_type), POINTER :: perturbed_density, rho
134 : TYPE(section_vals_type), POINTER :: input, xc_section
135 : TYPE(xc_derivative_set_type) :: deriv_set
136 : TYPE(xc_rho_set_type) :: rho_set
137 :
138 : ! Build the perturbed density matrix correction depending on the overlap derivative
139 : ! P1 = C0 C1 + C1 C0
140 : ! - C0_(mu j) S1_(jk) C0_(k nu)
141 : ! This routine is adapted from apply_op_2_dft. There, build_dm_response builds
142 : ! C0 * dCR + dCR * C0.
143 : ! build_dm_response is computing $-1 * (C^0 C^1 + C^1 C^0)$ and later on in the
144 : ! integration the factor 2 is applied to account for the occupancy.
145 : ! The sign is negative because the kernel is on the RHS of the Sternheimer equation.
146 : !
147 : ! The correction factor in this routine needs to have
148 : ! the opposite sign mathematically as (C0 C1 + C1 C0)
149 : ! so the same sign in the code because of the $-1$ in dCR
150 : ! so the opposite sign in the code because we are on the LHS of the Sternheimer equation.
151 : !
152 : ! This term must not go into the kernel applied by the linear response solver, because
153 : ! for the (P)CG algorithm, all constant terms have to be on one side of the equations
154 : ! and all solution dependent terms must be on the other side.
155 :
156 252 : CALL timeset(routineN, handle)
157 :
158 252 : NULLIFY (auxbas_pw_pool, pw_env, rho1_r, rho1_g_pw, &
159 252 : v_xc, poisson_env, input, rho, rho1_g, v_xc_tau)
160 :
161 252 : CALL cp_fm_create(rho_ao_fm, dcdr_env%aoao_fm_struct)
162 252 : CALL cp_fm_create(rho_ao_s1, dcdr_env%aoao_fm_struct)
163 252 : CALL cp_fm_create(rho_ao_s1_rho_ao, dcdr_env%aoao_fm_struct)
164 252 : CALL cp_fm_create(s1_ao, dcdr_env%aoao_fm_struct)
165 :
166 252 : IF (PRESENT(overlap1)) THEN
167 0 : CALL copy_dbcsr_to_fm(overlap1%matrix, s1_ao)
168 : ELSE
169 252 : CALL copy_dbcsr_to_fm(dcdr_env%matrix_s1(dcdr_env%beta + 1)%matrix, s1_ao)
170 : END IF
171 :
172 576 : DO ispin = 1, dcdr_env%nspins
173 324 : CALL dbcsr_set(dcdr_env%perturbed_dm_correction(ispin)%matrix, 0._dp)
174 324 : CALL dbcsr_set(dcdr_env%matrix_apply_op_constant(ispin)%matrix, 0.0_dp)
175 :
176 : CALL parallel_gemm('N', 'T', dcdr_env%nao, dcdr_env%nao, dcdr_env%nmo(ispin), &
177 : 1.0_dp, dcdr_env%mo_coeff(ispin), dcdr_env%mo_coeff(ispin), &
178 324 : 0.0_dp, rho_ao_fm)
179 :
180 : CALL parallel_gemm('N', 'N', dcdr_env%nao, dcdr_env%nao, dcdr_env%nao, &
181 : 1.0_dp, rho_ao_fm, s1_ao, &
182 324 : 0.0_dp, rho_ao_s1)
183 :
184 : CALL parallel_gemm('N', 'N', dcdr_env%nao, dcdr_env%nao, dcdr_env%nao, &
185 : -1._dp, rho_ao_s1, rho_ao_fm, & ! this is the sign mentioned above.
186 324 : 0.0_dp, rho_ao_s1_rho_ao)
187 :
188 576 : CALL copy_fm_to_dbcsr(rho_ao_s1_rho_ao, dcdr_env%perturbed_dm_correction(ispin)%matrix)
189 : END DO
190 :
191 252 : CALL cp_fm_release(rho_ao_fm)
192 252 : CALL cp_fm_release(rho_ao_s1)
193 252 : CALL cp_fm_release(rho_ao_s1_rho_ao)
194 252 : CALL cp_fm_release(s1_ao)
195 : ! Done building the density matrix correction
196 :
197 : ! Build the density struct from the environment
198 252 : NULLIFY (perturbed_density)
199 252 : ALLOCATE (perturbed_density)
200 252 : CALL qs_rho_create(perturbed_density)
201 252 : CALL qs_rho_rebuild(perturbed_density, qs_env=qs_env)
202 :
203 : ! ... set the density matrix to be the perturbed density matrix
204 252 : CALL qs_rho_get(perturbed_density, rho_ao=rho1_ao)
205 576 : DO ispin = 1, dcdr_env%nspins
206 576 : CALL dbcsr_copy(rho1_ao(ispin)%matrix, dcdr_env%perturbed_dm_correction(ispin)%matrix)
207 : END DO
208 :
209 : ! ... updates rho_r and rho_g to the rho%rho_ao.
210 : CALL qs_rho_update_rho(rho_struct=perturbed_density, &
211 252 : qs_env=qs_env)
212 :
213 : ! Also update the qs_env%rho
214 252 : CALL get_qs_env(qs_env, rho=rho)
215 252 : CALL qs_rho_update_rho(rho, qs_env=qs_env)
216 252 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
217 :
218 : energy_hartree = 0.0_dp
219 :
220 : CALL get_qs_env(qs_env=qs_env, &
221 : pw_env=pw_env, &
222 : xcint_weights=weights, &
223 252 : input=input)
224 :
225 : ! Create the temporary grids
226 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
227 252 : poisson_env=poisson_env)
228 :
229 : ! Allocate deriv_set and rho_set
230 252 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
231 :
232 : CALL xc_prep_2nd_deriv(deriv_set, rho_set, &
233 : rho_r, auxbas_pw_pool, weights, &
234 252 : xc_section=xc_section)
235 :
236 : ! Done with deriv_set and rho_set
237 :
238 1080 : ALLOCATE (v_rspace_new(dcdr_env%nspins))
239 252 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
240 252 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
241 :
242 : ! Calculate the Hartree potential on the total density
243 252 : CALL auxbas_pw_pool%create_pw(rho1_tot_gspace)
244 :
245 252 : CALL qs_rho_get(perturbed_density, rho_g=rho1_g, rho_r=rho1_r, tau_r=tau1_r)
246 252 : CALL pw_copy(rho1_g(1), rho1_tot_gspace)
247 324 : DO ispin = 2, dcdr_env%nspins
248 324 : CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
249 : END DO
250 :
251 : CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &
252 : energy_hartree, &
253 252 : v_hartree_gspace)
254 252 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
255 :
256 252 : CALL auxbas_pw_pool%give_back_pw(rho1_tot_gspace)
257 :
258 : ! Calculate the second derivative of the exchange-correlation potential
259 : CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, deriv_set, rho_set, &
260 : rho1_r, rho1_g_pw, tau1_r, auxbas_pw_pool, &
261 252 : weights, xc_section, gapw=.FALSE.)
262 :
263 576 : DO ispin = 1, dcdr_env%nspins
264 576 : v_rspace_new(ispin) = v_xc(ispin)
265 : END DO
266 252 : DEALLOCATE (v_xc)
267 :
268 : ! Done calculating the potentials
269 :
270 : !-------------------------------!
271 : ! Add both hartree and xc terms !
272 : !-------------------------------!
273 252 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
274 576 : DO ispin = 1, dcdr_env%nspins
275 576 : CALL pw_scale(v_rspace_new(ispin), v_rspace_new(ispin)%pw_grid%dvol)
276 : END DO
277 :
278 576 : DO ispin = 1, dcdr_env%nspins
279 324 : CALL dbcsr_set(dcdr_env%matrix_apply_op_constant(ispin)%matrix, 0.0_dp)
280 324 : CALL pw_axpy(v_hartree_rspace, v_rspace_new(ispin))
281 324 : IF (dcdr_env%nspins == 1) THEN
282 180 : CALL pw_scale(v_rspace_new(1), 2.0_dp)
283 : END IF
284 :
285 : CALL integrate_v_rspace(v_rspace=v_rspace_new(ispin), &
286 : hmat=dcdr_env%matrix_apply_op_constant(ispin), &
287 : qs_env=qs_env, &
288 576 : calculate_forces=.FALSE.)
289 : END DO
290 :
291 252 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
292 252 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
293 576 : DO ispin = 1, dcdr_env%nspins
294 576 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
295 : END DO
296 252 : DEALLOCATE (v_rspace_new)
297 :
298 252 : IF (ASSOCIATED(v_xc_tau)) THEN
299 0 : CALL pw_scale(v_xc_tau(1), 2._dp*v_xc_tau(1)%pw_grid%dvol)
300 : CALL integrate_v_rspace(v_rspace=v_xc_tau(1), &
301 : hmat=dcdr_env%matrix_apply_op_constant(1), &
302 : qs_env=qs_env, &
303 : compute_tau=.TRUE., &
304 0 : calculate_forces=.FALSE.)
305 :
306 0 : CALL auxbas_pw_pool%give_back_pw(v_xc_tau(1))
307 0 : DEALLOCATE (v_xc_tau)
308 : END IF
309 :
310 252 : CALL qs_rho_release(perturbed_density)
311 252 : DEALLOCATE (perturbed_density)
312 252 : CALL xc_rho_set_release(rho_set, auxbas_pw_pool)
313 252 : CALL xc_dset_release(deriv_set)
314 :
315 252 : CALL timestop(handle)
316 :
317 6048 : END SUBROUTINE apply_op_constant_term
318 :
319 : ! **************************************************************************************************
320 : !> \brief Calculate the derivative of the Hartree term due to the core charge density
321 : !> \param qs_env ...
322 : !> \param dcdr_env ...
323 : !> \author Edward Ditler
324 : ! **************************************************************************************************
325 72 : SUBROUTINE d_core_charge_density_dR(qs_env, dcdr_env)
326 : ! drho_core contribution
327 : ! sum over all directions
328 : ! output in ao x ao
329 : TYPE(qs_environment_type), POINTER :: qs_env
330 : TYPE(dcdr_env_type) :: dcdr_env
331 :
332 : CHARACTER(len=*), PARAMETER :: routineN = 'd_core_charge_density_dR'
333 :
334 : INTEGER :: beta, handle
335 : TYPE(cp_logger_type), POINTER :: logger
336 : TYPE(dft_control_type), POINTER :: dft_control
337 : TYPE(pw_c1d_gs_type) :: drho_g, v_hartree_gspace
338 : TYPE(pw_env_type), POINTER :: pw_env
339 : TYPE(pw_poisson_type), POINTER :: poisson_env
340 72 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
341 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
342 : TYPE(pw_r3d_rs_type) :: v_hartree_rspace
343 : TYPE(qs_rho_type), POINTER :: rho
344 :
345 72 : CALL timeset(routineN, handle)
346 :
347 72 : logger => cp_get_default_logger()
348 :
349 72 : NULLIFY (pw_env, auxbas_pw_pool, pw_pools, poisson_env, dft_control, &
350 72 : rho)
351 :
352 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho, &
353 72 : dft_control=dft_control)
354 :
355 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, poisson_env=poisson_env, &
356 72 : pw_pools=pw_pools)
357 :
358 : ! Create the Hartree potential grids in real and reciprocal space.
359 72 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
360 72 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
361 : ! Create the grid for the derivative of the core potential
362 72 : CALL auxbas_pw_pool%create_pw(drho_g)
363 :
364 288 : DO beta = 1, 3
365 216 : CALL pw_zero(v_hartree_gspace)
366 216 : CALL pw_zero(v_hartree_rspace)
367 216 : CALL pw_zero(drho_g)
368 :
369 : ! Calculate the Hartree potential on the perturbed density and Poisson solve it
370 : CALL calculate_drho_core(drho_core=drho_g, qs_env=qs_env, &
371 216 : beta=beta, lambda=dcdr_env%lambda)
372 : CALL pw_poisson_solve(poisson_env, drho_g, &
373 216 : vhartree=v_hartree_gspace)
374 216 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
375 216 : CALL pw_scale(v_hartree_rspace, v_hartree_rspace%pw_grid%dvol)
376 :
377 : ! Calculate the integrals
378 : CALL integrate_v_rspace(v_rspace=v_hartree_rspace, &
379 : hmat=dcdr_env%matrix_core_charge_1(beta), &
380 : qs_env=qs_env, &
381 288 : calculate_forces=.FALSE.)
382 : END DO
383 :
384 72 : CALL auxbas_pw_pool%give_back_pw(drho_g)
385 72 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
386 72 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
387 :
388 72 : CALL timestop(handle)
389 72 : END SUBROUTINE d_core_charge_density_dR
390 :
391 : ! **************************************************************************************************
392 : !> \brief Core Hamiltonian contributions to the operator (the pseudopotentials)
393 : !> \param qs_env ...
394 : !> \param dcdr_env ..
395 : !> \author Edward Ditler
396 : ! **************************************************************************************************
397 72 : SUBROUTINE core_dR(qs_env, dcdr_env)
398 : TYPE(qs_environment_type), POINTER :: qs_env
399 : TYPE(dcdr_env_type) :: dcdr_env
400 :
401 : CHARACTER(LEN=*), PARAMETER :: routineN = 'core_dR'
402 :
403 : CHARACTER(LEN=default_string_length) :: my_basis_type
404 : INTEGER :: handle, nder
405 : LOGICAL :: calculate_forces
406 72 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
407 72 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_h, matrix_p_pass
408 : TYPE(qs_ks_env_type), POINTER :: ks_env
409 : TYPE(qs_rho_type), POINTER :: rho
410 :
411 72 : CALL timeset(routineN, handle)
412 :
413 72 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env)
414 72 : CALL get_ks_env(ks_env=ks_env, rho=rho)
415 72 : CALL qs_rho_get(rho, rho_ao=rho_ao)
416 :
417 72 : nder = 1
418 72 : calculate_forces = .FALSE.
419 :
420 : my_basis_type = "ORB"
421 :
422 72 : NULLIFY (matrix_h)
423 72 : matrix_p_pass(1:1, 1:1) => rho_ao(1:1)
424 : CALL core_matrices(qs_env, matrix_h, matrix_p_pass, calculate_forces, nder, &
425 72 : dcdr_env=dcdr_env)
426 :
427 72 : CALL timestop(handle)
428 :
429 72 : END SUBROUTINE core_dR
430 :
431 : ! **************************************************************************************************
432 : !> \brief The derivatives of the basis functions going into the HXC potential wrt nuclear positions
433 : !> \param qs_env ...
434 : !> \param dcdr_env ...
435 : !> \author Edward Ditler
436 : ! **************************************************************************************************
437 72 : SUBROUTINE d_vhxc_dR(qs_env, dcdr_env)
438 : TYPE(qs_environment_type), POINTER :: qs_env
439 : TYPE(dcdr_env_type) :: dcdr_env
440 :
441 : CHARACTER(len=*), PARAMETER :: routineN = 'd_vhxc_dR'
442 :
443 : INTEGER :: handle, idir, ispin
444 72 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rho_ao
445 : TYPE(pw_c1d_gs_type) :: drho_g_total, v_hartree_gspace
446 72 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: drho_g
447 : TYPE(pw_env_type), POINTER :: pw_env
448 : TYPE(pw_poisson_type), POINTER :: poisson_env
449 72 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
450 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
451 : TYPE(pw_r3d_rs_type) :: drho_r_total, v_hartree_rspace
452 144 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: drho_r, dtau_r, rho_r, v_xc, v_xc_tau
453 : TYPE(pw_r3d_rs_type), POINTER :: weights
454 : TYPE(qs_rho_type), POINTER :: rho
455 : TYPE(section_vals_type), POINTER :: input, xc_section
456 : TYPE(xc_derivative_set_type) :: my_deriv_set
457 : TYPE(xc_rho_set_type) :: my_rho_set
458 :
459 72 : CALL timeset(routineN, handle)
460 :
461 : CALL get_qs_env(qs_env=qs_env, &
462 : pw_env=pw_env, &
463 : input=input, &
464 : xcint_weights=weights, &
465 72 : rho=rho)
466 72 : CALL qs_rho_get(rho, rho_ao=rho_ao, rho_r=rho_r)
467 :
468 72 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
469 :
470 : ! get the tmp grids
471 300 : ALLOCATE (drho_r(dcdr_env%nspins))
472 300 : ALLOCATE (drho_g(dcdr_env%nspins))
473 :
474 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
475 72 : pw_pools=pw_pools, poisson_env=poisson_env)
476 72 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
477 72 : CALL auxbas_pw_pool%create_pw(v_hartree_rspace)
478 :
479 156 : DO ispin = 1, dcdr_env%nspins
480 84 : CALL auxbas_pw_pool%create_pw(drho_r(ispin))
481 156 : CALL auxbas_pw_pool%create_pw(drho_g(ispin))
482 : END DO
483 72 : CALL auxbas_pw_pool%create_pw(drho_g_total)
484 72 : CALL auxbas_pw_pool%create_pw(drho_r_total)
485 :
486 288 : DO idir = 1, 3
487 216 : CALL pw_zero(v_hartree_gspace)
488 216 : CALL pw_zero(v_hartree_rspace)
489 216 : CALL pw_zero(drho_g_total)
490 216 : CALL pw_zero(drho_r_total)
491 :
492 468 : DO ispin = 1, dcdr_env%nspins
493 252 : CALL pw_zero(drho_r(ispin))
494 252 : CALL pw_zero(drho_g(ispin))
495 :
496 : ! Get the density
497 : CALL calculate_drho_elec_dR(matrix_p=rho_ao(ispin)%matrix, &
498 : drho=drho_r(ispin), &
499 : drho_gspace=drho_g(ispin), &
500 : qs_env=qs_env, &
501 252 : beta=idir, lambda=dcdr_env%lambda)
502 :
503 252 : CALL pw_axpy(drho_g(ispin), drho_g_total)
504 468 : CALL pw_axpy(drho_r(ispin), drho_r_total)
505 : END DO
506 : ! Get the Hartree potential corresponding to the perturbed density
507 : CALL pw_poisson_solve(poisson_env, drho_g_total, &
508 216 : vhartree=v_hartree_gspace)
509 216 : CALL pw_transfer(v_hartree_gspace, v_hartree_rspace)
510 :
511 : ! Get the XC potential corresponding to the perturbed density
512 : CALL xc_prep_2nd_deriv(my_deriv_set, my_rho_set, &
513 : rho_r, auxbas_pw_pool, weights, &
514 216 : xc_section=xc_section)
515 :
516 216 : NULLIFY (dtau_r)
517 : CALL xc_calc_2nd_deriv(v_xc, v_xc_tau, my_deriv_set, my_rho_set, &
518 : drho_r, drho_g, dtau_r, auxbas_pw_pool, &
519 216 : weights, xc_section, gapw=.FALSE.)
520 216 : IF (ASSOCIATED(v_xc_tau)) THEN
521 0 : CPABORT(dcdr_meta_gga_error)
522 : END IF
523 :
524 216 : CALL xc_dset_release(my_deriv_set)
525 216 : CALL xc_rho_set_release(my_rho_set)
526 :
527 : !-------------------------------!
528 : ! Add both hartree and xc terms !
529 : !-------------------------------!
530 468 : DO ispin = 1, dcdr_env%nspins
531 : ! Can the dvol be different?
532 252 : CALL pw_scale(v_xc(ispin), v_xc(ispin)%pw_grid%dvol)
533 252 : CALL pw_axpy(v_hartree_rspace, v_xc(ispin), v_hartree_rspace%pw_grid%dvol)
534 :
535 : CALL integrate_v_rspace(v_rspace=v_xc(ispin), &
536 : hmat=dcdr_env%matrix_d_vhxc_dR(idir, ispin), &
537 : qs_env=qs_env, &
538 252 : calculate_forces=.FALSE.)
539 :
540 : ! v_xc gets allocated again in xc_calc_2nd_deriv
541 468 : CALL auxbas_pw_pool%give_back_pw(v_xc(ispin))
542 : END DO ! ispin
543 288 : DEALLOCATE (v_xc)
544 : END DO ! idir
545 :
546 72 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
547 72 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rspace)
548 72 : CALL auxbas_pw_pool%give_back_pw(drho_g_total)
549 72 : CALL auxbas_pw_pool%give_back_pw(drho_r_total)
550 :
551 156 : DO ispin = 1, dcdr_env%nspins
552 84 : CALL auxbas_pw_pool%give_back_pw(drho_g(ispin))
553 156 : CALL auxbas_pw_pool%give_back_pw(drho_r(ispin))
554 : END DO
555 :
556 72 : DEALLOCATE (drho_g)
557 72 : DEALLOCATE (drho_r)
558 :
559 72 : CALL timestop(handle)
560 :
561 1584 : END SUBROUTINE d_vhxc_dR
562 :
563 : ! **************************************************************************************************
564 : !> \brief The derivatives of the basis functions over which the HXC potential is integrated,
565 : !> so < da/dR | Vhxc | b >
566 : !> \param qs_env ...
567 : !> \param dcdr_env ...
568 : !> \author Edward Ditler
569 : ! **************************************************************************************************
570 72 : SUBROUTINE vhxc_R_perturbed_basis_functions(qs_env, dcdr_env)
571 : TYPE(qs_environment_type), POINTER :: qs_env
572 : TYPE(dcdr_env_type) :: dcdr_env
573 :
574 : CHARACTER(LEN=*), PARAMETER :: routineN = 'vhxc_R_perturbed_basis_functions'
575 :
576 : INTEGER :: handle, ispin
577 72 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_vhxc_dbasis
578 72 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
579 : TYPE(pw_env_type), POINTER :: pw_env
580 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
581 72 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: v_hxc_r, v_tau_rspace
582 : TYPE(pw_r3d_rs_type), POINTER :: v_hartree_r
583 : TYPE(qs_energy_type), POINTER :: energy
584 : TYPE(qs_ks_env_type), POINTER :: ks_env
585 : TYPE(qs_rho_type), POINTER :: rho_struct
586 : TYPE(section_vals_type), POINTER :: input, xc_section
587 :
588 72 : CALL timeset(routineN, handle)
589 :
590 72 : NULLIFY (rho_struct, energy, input, ks_env, pw_env, matrix_p)
591 : CALL get_qs_env(qs_env, &
592 : rho=rho_struct, &
593 : energy=energy, &
594 : input=input, &
595 : ks_env=ks_env, &
596 : pw_env=pw_env, &
597 72 : v_hartree_rspace=v_hartree_r)
598 72 : CALL qs_rho_get(rho_struct, rho_ao_kp=matrix_p)
599 72 : xc_section => section_vals_get_subs_vals(input, "DFT%XC")
600 :
601 72 : NULLIFY (auxbas_pw_pool)
602 72 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
603 :
604 : ! *** calculate the xc potential on the pw density ***
605 : ! *** associates v_hxc_r if the xc potential needs to be computed.
606 : ! If we do wavefunction fitting, we need the vxc_potential in the auxiliary basis set
607 72 : NULLIFY (v_hxc_r, v_tau_rspace)
608 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_struct, xc_section=xc_section, &
609 72 : vxc_rho=v_hxc_r, vxc_tau=v_tau_rspace, exc=energy%exc)
610 :
611 156 : DO ispin = 1, dcdr_env%nspins
612 84 : CALL pw_scale(v_hxc_r(ispin), v_hxc_r(ispin)%pw_grid%dvol)
613 :
614 : ! sum up potentials and integrate
615 84 : CALL pw_axpy(v_hartree_r, v_hxc_r(ispin), 1._dp)
616 :
617 84 : matrix_vhxc_dbasis => dcdr_env%matrix_vhxc_perturbed_basis(ispin, :)
618 : CALL integrate_v_dbasis(v_rspace=v_hxc_r(ispin), &
619 : matrix_p=matrix_p(ispin, 1)%matrix, &
620 : matrix_vhxc_dbasis=matrix_vhxc_dbasis, &
621 : qs_env=qs_env, &
622 84 : lambda=dcdr_env%lambda)
623 :
624 156 : CALL auxbas_pw_pool%give_back_pw(v_hxc_r(ispin))
625 : END DO
626 :
627 72 : DEALLOCATE (v_hxc_r)
628 :
629 72 : CALL timestop(handle)
630 72 : END SUBROUTINE vhxc_R_perturbed_basis_functions
631 :
632 : ! **************************************************************************************************
633 : !> \brief Enforce that one of the basis functions in < a | O | b > is centered on atom lambda.
634 : !> \param matrix ...
635 : !> \param qs_kind_set ...
636 : !> \param basis_type ...
637 : !> \param sab_nl ...
638 : !> \param lambda Atom index
639 : !> \param direction_Or True: < a | O | b==lambda >, False: < a==lambda | O | b >
640 : ! **************************************************************************************************
641 2610 : SUBROUTINE hr_mult_by_delta_1d(matrix, qs_kind_set, basis_type, sab_nl, lambda, direction_Or)
642 :
643 : TYPE(dbcsr_type), POINTER :: matrix
644 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
645 : CHARACTER(LEN=*), INTENT(IN) :: basis_type
646 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
647 : POINTER :: sab_nl
648 : INTEGER, INTENT(IN) :: lambda
649 : LOGICAL, INTENT(IN) :: direction_Or
650 :
651 : CHARACTER(len=*), PARAMETER :: routineN = 'hr_mult_by_delta_1d'
652 :
653 : INTEGER :: handle, iatom, icol, ikind, irow, jatom, &
654 : jkind, ldsab, mepos, nkind, nseta, &
655 : nsetb, nthread
656 : INTEGER, DIMENSION(3) :: cell
657 2610 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
658 2610 : npgfb, nsgfa, nsgfb
659 2610 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
660 : LOGICAL :: do_symmetric, found
661 : REAL(KIND=dp), DIMENSION(3) :: rab
662 2610 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a, set_radius_b
663 2610 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: k_block, rpgfa, rpgfb, scon_a, scon_b, &
664 2610 : zeta, zetb
665 2610 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
666 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
667 : TYPE(neighbor_list_iterator_p_type), &
668 2610 : DIMENSION(:), POINTER :: nl_iterator
669 :
670 2610 : CALL timeset(routineN, handle)
671 :
672 2610 : nkind = SIZE(qs_kind_set)
673 :
674 : ! check for symmetry
675 2610 : CPASSERT(SIZE(sab_nl) > 0)
676 2610 : CALL get_neighbor_list_set_p(neighbor_list_sets=sab_nl, symmetric=do_symmetric)
677 :
678 : ! prepare basis set
679 13050 : ALLOCATE (basis_set_list(nkind))
680 2610 : CALL basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
681 :
682 : ! *** Allocate work storage ***
683 2610 : ldsab = get_memory_usage(qs_kind_set, basis_type)
684 :
685 : nthread = 1
686 2610 : !$ nthread = omp_get_max_threads()
687 : ! Iterate of neighbor list
688 2610 : CALL neighbor_list_iterator_create(nl_iterator, sab_nl, nthread=nthread)
689 :
690 : !$OMP PARALLEL DEFAULT(NONE) &
691 : !$OMP SHARED (nthread,ldsab,nl_iterator, do_symmetric) &
692 : !$OMP SHARED (ncoset,matrix,basis_set_list) &
693 : !$OMP SHARED (direction_or, lambda) &
694 : !$OMP PRIVATE (k_block,mepos,ikind,jkind,iatom,jatom,rab,cell) &
695 : !$OMP PRIVATE (basis_set_a,basis_set_b) &
696 : !$OMP PRIVATE (first_sgfa, la_max, la_min, npgfa, nsgfa, nseta, rpgfa, set_radius_a) &
697 : !$OMP PRIVATE (zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, nsgfb) &
698 2610 : !$OMP PRIVATE (zetb, scon_a, scon_b, irow, icol, found)
699 :
700 : mepos = 0
701 : !$ mepos = omp_get_thread_num()
702 :
703 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
704 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, &
705 : iatom=iatom, jatom=jatom, r=rab, cell=cell)
706 : basis_set_a => basis_set_list(ikind)%gto_basis_set
707 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
708 : basis_set_b => basis_set_list(jkind)%gto_basis_set
709 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
710 : ! basis ikind
711 : first_sgfa => basis_set_a%first_sgf
712 : la_max => basis_set_a%lmax
713 : la_min => basis_set_a%lmin
714 : npgfa => basis_set_a%npgf
715 : nseta = basis_set_a%nset
716 : nsgfa => basis_set_a%nsgf_set
717 : rpgfa => basis_set_a%pgf_radius
718 : set_radius_a => basis_set_a%set_radius
719 : scon_a => basis_set_a%scon
720 : zeta => basis_set_a%zet
721 : ! basis jkind
722 : first_sgfb => basis_set_b%first_sgf
723 : lb_max => basis_set_b%lmax
724 : lb_min => basis_set_b%lmin
725 : npgfb => basis_set_b%npgf
726 : nsetb = basis_set_b%nset
727 : nsgfb => basis_set_b%nsgf_set
728 : rpgfb => basis_set_b%pgf_radius
729 : set_radius_b => basis_set_b%set_radius
730 : scon_b => basis_set_b%scon
731 : zetb => basis_set_b%zet
732 :
733 : IF (do_symmetric) THEN
734 : IF (iatom <= jatom) THEN
735 : irow = iatom
736 : icol = jatom
737 : ELSE
738 : irow = jatom
739 : icol = iatom
740 : END IF
741 : ELSE
742 : irow = iatom
743 : icol = jatom
744 : END IF
745 :
746 : NULLIFY (k_block)
747 : CALL dbcsr_get_block_p(matrix, irow, icol, k_block, found)
748 : CPASSERT(found)
749 :
750 : IF (direction_Or) THEN
751 : IF (jatom /= lambda) k_block(:, :) = 0._dp
752 : ELSE IF (.NOT. direction_Or) THEN
753 : IF (iatom /= lambda) k_block(:, :) = 0._dp
754 : END IF
755 : END DO
756 : !$OMP END PARALLEL
757 2610 : CALL neighbor_list_iterator_release(nl_iterator)
758 :
759 : ! Release work storage
760 2610 : DEALLOCATE (basis_set_list)
761 :
762 2610 : CALL timestop(handle)
763 :
764 5220 : END SUBROUTINE hr_mult_by_delta_1d
765 :
766 : END MODULE qs_dcdr_ao
|
