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