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