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 Harris method calculations
10 : !> \par History
11 : !> 2026.07 split from qs_harris_utils
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_harris_methods
15 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
16 : dbcsr_create,&
17 : dbcsr_p_type,&
18 : dbcsr_release,&
19 : dbcsr_set
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_get_default_io_unit,&
22 : cp_logger_type
23 : USE distribution_1d_types, ONLY: distribution_1d_type
24 : USE input_constants, ONLY: hden_cube,&
25 : hden_cube_fit
26 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
27 : section_vals_type
28 : USE kinds, ONLY: dp
29 : USE message_passing, ONLY: mp_para_env_type
30 : USE pw_env_types, ONLY: pw_env_get,&
31 : pw_env_type
32 : USE pw_methods, ONLY: pw_axpy,&
33 : pw_copy,&
34 : pw_integral_ab,&
35 : pw_integrate_function,&
36 : pw_scale,&
37 : pw_transfer
38 : USE pw_poisson_methods, ONLY: pw_poisson_solve
39 : USE pw_poisson_types, ONLY: pw_poisson_type
40 : USE pw_pool_types, ONLY: pw_pool_type
41 : USE pw_types, ONLY: pw_c1d_gs_type,&
42 : pw_r3d_rs_type
43 : USE qs_collocate_density, ONLY: calculate_rho_elec
44 : USE qs_density_fit, ONLY: fit_relative_entropy_density
45 : USE qs_energy_types, ONLY: qs_energy_type
46 : USE qs_environment_types, ONLY: get_qs_env,&
47 : qs_environment_type
48 : USE qs_force_types, ONLY: qs_force_type
49 : USE qs_fxc, ONLY: qs_fxc_create
50 : USE qs_harris_types, ONLY: harris_energy_type,&
51 : harris_print_direct_energy,&
52 : harris_print_energy,&
53 : harris_rhoin_type,&
54 : harris_type
55 : USE qs_integrate_potential, ONLY: integrate_function,&
56 : integrate_v_core_rspace,&
57 : integrate_v_rspace
58 : USE qs_ks_methods, ONLY: qs_ks_update_qs_env
59 : USE qs_ks_types, ONLY: qs_ks_did_change,&
60 : qs_ks_env_type
61 : USE qs_rho_atom_types, ONLY: rho_atom_type
62 : USE qs_rho_types, ONLY: qs_rho_create,&
63 : qs_rho_get,&
64 : qs_rho_set,&
65 : qs_rho_type
66 : #include "./base/base_uses.f90"
67 :
68 : IMPLICIT NONE
69 :
70 : PRIVATE
71 :
72 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_harris_methods'
73 :
74 : PUBLIC :: harris_direct_density_matrix_energy, harris_energy_correction, &
75 : harris_relative_entropy_reconstruction
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief ...
81 : !> \param qs_env ...
82 : !> \param calculate_forces ...
83 : ! **************************************************************************************************
84 64 : SUBROUTINE harris_energy_correction(qs_env, calculate_forces)
85 : TYPE(qs_environment_type), POINTER :: qs_env
86 : LOGICAL, INTENT(IN) :: calculate_forces
87 :
88 : CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_energy_correction'
89 :
90 : INTEGER :: handle, iounit, ispin, nspins
91 : REAL(KIND=dp) :: dvol, ec, eh, exc, vxc
92 : TYPE(cp_logger_type), POINTER :: logger
93 : TYPE(harris_energy_type), POINTER :: energy
94 : TYPE(harris_type), POINTER :: harris_env
95 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
96 : TYPE(pw_env_type), POINTER :: pw_env
97 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
98 : TYPE(pw_r3d_rs_type) :: core_rspace
99 64 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
100 : TYPE(qs_energy_type), POINTER :: ks_energy
101 : TYPE(qs_rho_type), POINTER :: rho
102 :
103 : MARK_USED(calculate_forces)
104 :
105 64 : CALL timeset(routineN, handle)
106 :
107 64 : CALL get_qs_env(qs_env, harris_env=harris_env, energy=ks_energy)
108 64 : energy => harris_env%energy
109 64 : IF (harris_env%direct_density_matrix_energy) THEN
110 8 : IF (calculate_forces) THEN
111 0 : CPABORT("Forces are not available for the direct fitted-density-matrix energy")
112 : END IF
113 8 : ks_energy%nonscf_correction = energy%direct_harris - ks_energy%total
114 8 : ks_energy%total = energy%direct_harris
115 8 : CALL timestop(handle)
116 8 : RETURN
117 : END IF
118 56 : IF (calculate_forces .AND. (harris_env%density_source == hden_cube .OR. &
119 : harris_env%density_source == hden_cube_fit)) THEN
120 0 : CPABORT("Forces are not available for a Harris energy based on an external cube density")
121 : END IF
122 56 : energy%eband = ks_energy%band
123 56 : energy%ewald_correction = ks_energy%core_overlap + ks_energy%core_self
124 56 : energy%dispersion = ks_energy%dispersion
125 :
126 56 : nspins = harris_env%rhoin%nspin
127 :
128 56 : CALL get_qs_env(qs_env, rho=rho, rho_core=rho_core)
129 56 : CALL qs_rho_get(rho, rho_r=rho_r)
130 :
131 56 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
132 56 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
133 56 : CALL auxbas_pw_pool%create_pw(core_rspace)
134 56 : CALL pw_transfer(rho_core, core_rspace)
135 :
136 56 : dvol = harris_env%vh_rspace%pw_grid%dvol
137 56 : eh = 0.0_dp
138 126 : DO ispin = 1, nspins
139 126 : eh = eh + pw_integral_ab(rho_r(ispin), harris_env%vh_rspace)/dvol
140 : END DO
141 56 : ec = pw_integral_ab(core_rspace, harris_env%vh_rspace)/dvol
142 56 : eh = 0.5_dp*(eh + ec)
143 56 : energy%eh_correction = ec - eh
144 :
145 56 : exc = ks_energy%exc
146 56 : vxc = 0.0_dp
147 56 : IF (ASSOCIATED(harris_env%vxc_rspace)) THEN
148 126 : DO ispin = 1, nspins
149 : vxc = vxc + pw_integral_ab(rho_r(ispin), harris_env%vxc_rspace(ispin))/ &
150 126 : harris_env%vxc_rspace(ispin)%pw_grid%dvol
151 : END DO
152 : END IF
153 56 : energy%exc_correction = exc - vxc
154 :
155 : ! Total Harris model energy
156 : energy%eharris = energy%eband + energy%eh_correction + energy%exc_correction + &
157 56 : energy%ewald_correction + energy%dispersion
158 :
159 56 : CALL auxbas_pw_pool%give_back_pw(core_rspace)
160 :
161 56 : ks_energy%total = ks_energy%total + ks_energy%core
162 56 : ks_energy%nonscf_correction = energy%eharris - ks_energy%total
163 56 : ks_energy%total = energy%eharris
164 :
165 56 : logger => cp_get_default_logger()
166 56 : iounit = cp_logger_get_default_io_unit(logger)
167 :
168 56 : CALL harris_print_energy(iounit, energy)
169 :
170 56 : IF (calculate_forces) THEN
171 6 : CALL harris_forces(qs_env, iounit)
172 : END IF
173 :
174 56 : CALL timestop(handle)
175 :
176 64 : END SUBROUTINE harris_energy_correction
177 :
178 : ! **************************************************************************************************
179 : !> \brief Evaluates the fitted AO density matrix without solving a NONSCF eigenproblem.
180 : !>
181 : !> First, the fitted density matrix and its collocated density are used consistently. The second
182 : !> energy keeps all density-matrix-dependent terms from this evaluation and replaces only the
183 : !> Hartree and semilocal XC contributions by those obtained from the original cube density. For
184 : !> GPW this is algebraically equivalent to the Harris expression with P_fit in the band trace.
185 : !> \param qs_env Quickstep environment
186 : ! **************************************************************************************************
187 8 : SUBROUTINE harris_direct_density_matrix_energy(qs_env)
188 : TYPE(qs_environment_type), POINTER :: qs_env
189 :
190 : CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_direct_density_matrix_energy'
191 :
192 : INTEGER :: handle, iounit
193 : REAL(KIND=dp) :: target_density_energy, &
194 : trial_density_energy
195 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
196 : TYPE(cp_logger_type), POINTER :: logger
197 : TYPE(harris_energy_type), POINTER :: energy
198 : TYPE(harris_type), POINTER :: harris_env
199 8 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
200 8 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
201 : TYPE(qs_energy_type), POINTER :: ks_energy
202 : TYPE(qs_rho_type), POINTER :: rho
203 :
204 8 : CALL timeset(routineN, handle)
205 8 : NULLIFY (harris_env, ks_energy, rho, rho_g, rho_r, tot_rho_r)
206 :
207 8 : CALL get_qs_env(qs_env, harris_env=harris_env, energy=ks_energy, rho=rho)
208 8 : CPASSERT(harris_env%density_source == hden_cube_fit)
209 8 : CPASSERT(harris_env%density_fit_ready)
210 8 : CPASSERT(ASSOCIATED(harris_env%density_target_rspace%pw_grid))
211 :
212 8 : energy => harris_env%energy
213 :
214 : ! Consistent evaluation: P_fit and the density collocated from P_fit.
215 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE., &
216 8 : print_active=.FALSE.)
217 8 : energy%trial_dm = ks_energy%total
218 : trial_density_energy = ks_energy%hartree + ks_energy%hartree_1c + &
219 : ks_energy%exc + ks_energy%exc1 + &
220 8 : ks_energy%exc_aux_fit + ks_energy%exc1_aux_fit
221 :
222 : ! Mixed Harris evaluation: retain P_fit but use the original cube density for E_H and E_xc.
223 8 : CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
224 8 : CPASSERT(SIZE(rho_r) == 1 .AND. SIZE(rho_g) == 1)
225 8 : CALL pw_copy(harris_env%density_target_rspace, rho_r(1))
226 8 : CALL pw_transfer(rho_r(1), rho_g(1))
227 8 : tot_rho_r(1) = pw_integrate_function(rho_r(1), isign=-1)
228 8 : CALL qs_rho_set(rho, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
229 8 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE., potential_changed=.TRUE.)
230 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE., &
231 8 : print_active=.FALSE.)
232 : target_density_energy = ks_energy%hartree + ks_energy%hartree_1c + &
233 : ks_energy%exc + ks_energy%exc1 + &
234 8 : ks_energy%exc_aux_fit + ks_energy%exc1_aux_fit
235 :
236 8 : energy%direct_harris = energy%trial_dm + target_density_energy - trial_density_energy
237 8 : energy%direct_difference = energy%trial_dm - energy%direct_harris
238 8 : energy%eharris = energy%direct_harris
239 8 : ks_energy%nonscf_correction = energy%direct_harris - ks_energy%total
240 8 : ks_energy%total = energy%direct_harris
241 :
242 8 : logger => cp_get_default_logger()
243 8 : iounit = cp_logger_get_default_io_unit(logger)
244 8 : CALL harris_print_direct_energy(iounit, energy)
245 :
246 8 : CALL timestop(handle)
247 :
248 8 : END SUBROUTINE harris_direct_density_matrix_energy
249 :
250 : ! **************************************************************************************************
251 : !> \brief Builds H[n_cube] once and reconstructs the AO density matrix by minimizing a
252 : !> density residual regularized by fermionic relative entropy to its Fermi matrix.
253 : !> \param qs_env Quickstep environment with the target cube density active
254 : ! **************************************************************************************************
255 8 : SUBROUTINE harris_relative_entropy_reconstruction(qs_env)
256 : TYPE(qs_environment_type), POINTER :: qs_env
257 :
258 : CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_relative_entropy_reconstruction'
259 :
260 : INTEGER :: handle
261 8 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
262 8 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks
263 : TYPE(harris_type), POINTER :: harris_env
264 8 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
265 8 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
266 : TYPE(qs_rho_type), POINTER :: rho
267 :
268 8 : CALL timeset(routineN, handle)
269 8 : NULLIFY (harris_env, matrix_ks, rho, rho_g, rho_r, tot_rho_r)
270 :
271 8 : CALL get_qs_env(qs_env, harris_env=harris_env, rho=rho)
272 8 : CPASSERT(harris_env%density_source == hden_cube_fit)
273 8 : CPASSERT(harris_env%density_target_ready)
274 8 : CPASSERT(.NOT. harris_env%density_fit_ready)
275 :
276 : ! The active real-space density is n_cube here. Build and freeze its Kohn-Sham matrix
277 : ! before changing the AO density matrix.
278 8 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE., potential_changed=.TRUE.)
279 : CALL qs_ks_update_qs_env(qs_env, calculate_forces=.FALSE., just_energy=.TRUE., &
280 8 : print_active=.FALSE.)
281 8 : CALL get_qs_env(qs_env, matrix_ks=matrix_ks)
282 8 : CPASSERT(ASSOCIATED(matrix_ks) .AND. SIZE(matrix_ks) == 1)
283 :
284 : CALL fit_relative_entropy_density(qs_env, rho, matrix_ks(1)%matrix, &
285 : harris_env%fit_temperature, &
286 : harris_env%fit_relative_entropy_weight, &
287 : harris_env%fit_max_iter, harris_env%fit_eps, &
288 8 : harris_env%fit_step_size, harris_env%fit_max_backtrack)
289 :
290 8 : CALL qs_rho_get(rho, rho_r=rho_r, rho_g=rho_g, tot_rho_r=tot_rho_r)
291 8 : CPASSERT(SIZE(rho_r) == 1 .AND. SIZE(rho_g) == 1)
292 8 : CALL harris_env%density_fit_rspace%create(rho_r(1)%pw_grid)
293 8 : CALL pw_copy(rho_r(1), harris_env%density_fit_rspace)
294 8 : harris_env%density_fit_ready = .TRUE.
295 8 : CALL qs_rho_set(rho, rho_r_valid=.TRUE., rho_g_valid=.TRUE.)
296 8 : CALL qs_ks_did_change(qs_env%ks_env, rho_changed=.TRUE., potential_changed=.TRUE.)
297 :
298 8 : CALL timestop(handle)
299 :
300 8 : END SUBROUTINE harris_relative_entropy_reconstruction
301 :
302 : ! **************************************************************************************************
303 : !> \brief ...
304 : !> \param qs_env ...
305 : !> \param iounit ...
306 : ! **************************************************************************************************
307 6 : SUBROUTINE harris_forces(qs_env, iounit)
308 : TYPE(qs_environment_type), POINTER :: qs_env
309 : INTEGER, INTENT(IN) :: iounit
310 :
311 : CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_forces'
312 : LOGICAL, PARAMETER :: debug_forces = .TRUE.
313 :
314 : INTEGER :: handle, ispin, nspins
315 : REAL(KIND=dp) :: ehartree
316 : REAL(KIND=dp), DIMENSION(3) :: fodeb
317 : TYPE(dbcsr_p_type) :: scrm
318 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rhoh_ao, smat
319 : TYPE(harris_type), POINTER :: harris_env
320 : TYPE(mp_para_env_type), POINTER :: para_env
321 : TYPE(pw_c1d_gs_type) :: rhoh_tot_gspace, vhout_gspace
322 6 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rhoh_g
323 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
324 : TYPE(pw_env_type), POINTER :: pw_env
325 : TYPE(pw_poisson_type), POINTER :: poisson_env
326 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
327 : TYPE(pw_r3d_rs_type) :: vhout_rspace, vhxc_rspace
328 6 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fhxc_rspace, ftau, fxc, rho_r, rhoh_r, &
329 6 : tauh_r
330 6 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
331 : TYPE(qs_ks_env_type), POINTER :: ks_env
332 : TYPE(qs_rho_type), POINTER :: rho, rhoh
333 6 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set
334 : TYPE(section_vals_type), POINTER :: xc_section
335 :
336 6 : CALL timeset(routineN, handle)
337 :
338 : IF (debug_forces) THEN
339 6 : IF (iounit > 0) WRITE (iounit, "(/,T3,A)") &
340 3 : "DEBUG:: Harris Method Forces (density dependent)"
341 : END IF
342 :
343 6 : CALL get_qs_env(qs_env, harris_env=harris_env, force=force, para_env=para_env)
344 6 : nspins = harris_env%rhoin%nspin
345 :
346 6 : CALL get_qs_env(qs_env, rho=rho, rho_core=rho_core, matrix_s=smat)
347 : ! Warning: rho_ao = output DM; rho_r = rhoin
348 6 : CALL qs_rho_get(rho, rho_ao=rhoh_ao, rho_r=rho_r, rho_g=rho_g)
349 6 : ALLOCATE (scrm%matrix)
350 6 : CALL dbcsr_create(scrm%matrix, template=rhoh_ao(1)%matrix)
351 6 : CALL dbcsr_copy(scrm%matrix, smat(1)%matrix)
352 6 : CALL dbcsr_set(scrm%matrix, 0.0_dp)
353 :
354 6 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, ks_env=ks_env)
355 6 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
356 6 : CALL auxbas_pw_pool%create_pw(vhxc_rspace)
357 :
358 24 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
359 14 : DO ispin = 1, nspins
360 8 : CALL pw_copy(harris_env%vh_rspace, vhxc_rspace)
361 8 : CALL pw_axpy(harris_env%vxc_rspace(ispin), vhxc_rspace)
362 : CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
363 : hmat=scrm, pmat=rhoh_ao(ispin), &
364 14 : qs_env=qs_env, calculate_forces=.TRUE.)
365 : END DO
366 : IF (debug_forces) THEN
367 24 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
368 6 : CALL para_env%sum(fodeb)
369 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*(Vh[in]+Vxc)", fodeb
370 : END IF
371 :
372 6 : CALL dbcsr_release(scrm%matrix)
373 6 : DEALLOCATE (scrm%matrix)
374 6 : CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
375 :
376 46 : ALLOCATE (rhoh_r(nspins), rhoh_g(nspins))
377 14 : DO ispin = 1, nspins
378 8 : CALL auxbas_pw_pool%create_pw(rhoh_r(ispin))
379 14 : CALL auxbas_pw_pool%create_pw(rhoh_g(ispin))
380 : END DO
381 6 : CALL auxbas_pw_pool%create_pw(rhoh_tot_gspace)
382 6 : CALL pw_copy(rho_core, rhoh_tot_gspace)
383 14 : DO ispin = 1, nspins
384 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=rhoh_ao(ispin)%matrix, &
385 8 : rho=rhoh_r(ispin), rho_gspace=rhoh_g(ispin))
386 14 : CALL pw_axpy(rhoh_g(ispin), rhoh_tot_gspace)
387 : END DO
388 : ! no meta functionals here
389 6 : NULLIFY (tauh_r)
390 :
391 6 : CALL auxbas_pw_pool%create_pw(vhout_rspace)
392 6 : CALL auxbas_pw_pool%create_pw(vhout_gspace)
393 6 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
394 : !
395 6 : CALL pw_poisson_solve(poisson_env, rhoh_tot_gspace, ehartree, vhout_gspace)
396 : !
397 6 : CALL pw_transfer(vhout_gspace, vhout_rspace)
398 6 : CALL pw_scale(vhout_rspace, vhout_rspace%pw_grid%dvol)
399 :
400 24 : IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
401 6 : CALL integrate_v_core_rspace(vhout_rspace, qs_env)
402 : IF (debug_forces) THEN
403 24 : fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
404 6 : CALL para_env%sum(fodeb)
405 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh[out]*dncore ", fodeb
406 : END IF
407 :
408 20 : ALLOCATE (fhxc_rspace(nspins))
409 14 : DO ispin = 1, nspins
410 14 : CALL auxbas_pw_pool%create_pw(fhxc_rspace(ispin))
411 : END DO
412 : ! vh = vh[out] - vh[in]
413 6 : CALL pw_axpy(harris_env%vh_rspace, vhout_rspace, alpha=-1._dp, beta=1.0_dp)
414 : ! kernel fxc
415 : ! drho = rho[out] - rho[in]
416 14 : DO ispin = 1, nspins
417 8 : CALL pw_axpy(rho_r(ispin), rhoh_r(ispin), alpha=-1._dp, beta=1.0_dp)
418 14 : CALL pw_axpy(rho_g(ispin), rhoh_g(ispin), alpha=-1._dp, beta=1.0_dp)
419 : END DO
420 6 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
421 6 : NULLIFY (fxc, ftau)
422 6 : NULLIFY (rho0_atom_set, rho1_atom_set)
423 6 : ALLOCATE (rhoh)
424 6 : CALL qs_rho_create(rhoh)
425 6 : IF (ASSOCIATED(rhoh_r)) THEN
426 6 : CALL qs_rho_set(rhoh, rho_r=rhoh_r, rho_r_valid=.TRUE.)
427 : END IF
428 6 : IF (ASSOCIATED(rhoh_g)) THEN
429 6 : CALL qs_rho_set(rhoh, rho_g=rhoh_g, rho_g_valid=.TRUE.)
430 : END IF
431 6 : IF (ASSOCIATED(tauh_r)) THEN
432 0 : CALL qs_rho_set(rhoh, tau_r=tauh_r, tau_r_valid=.TRUE.)
433 : END IF
434 : !
435 : CALL qs_fxc_create(qs_env, rho, rhoh, rho0_atom_set, xc_section, .FALSE., &
436 6 : fxc, ftau, rho1_atom_set)
437 : !
438 6 : DEALLOCATE (rhoh)
439 6 : CPASSERT(.NOT. ASSOCIATED(ftau))
440 :
441 14 : DO ispin = 1, nspins
442 8 : CALL pw_copy(vhout_rspace, fhxc_rspace(ispin))
443 14 : IF (ASSOCIATED(fxc)) THEN
444 8 : CALL pw_scale(fxc(ispin), fxc(ispin)%pw_grid%dvol)
445 8 : CALL pw_axpy(fxc(ispin), fhxc_rspace(ispin))
446 : END IF
447 : END DO
448 :
449 24 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
450 6 : CALL calculate_harris_integrals(qs_env, harris_env%rhoin, fhxc_rspace, .TRUE.)
451 : IF (debug_forces) THEN
452 24 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
453 6 : CALL para_env%sum(fodeb)
454 6 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (dVh+fxc)*dn[in] ", fodeb
455 : END IF
456 :
457 6 : IF (ASSOCIATED(fxc)) THEN
458 14 : DO ispin = 1, nspins
459 14 : CALL auxbas_pw_pool%give_back_pw(fxc(ispin))
460 : END DO
461 6 : DEALLOCATE (fxc)
462 : END IF
463 6 : IF (ASSOCIATED(ftau)) THEN
464 0 : DO ispin = 1, nspins
465 0 : CALL auxbas_pw_pool%give_back_pw(ftau(ispin))
466 : END DO
467 0 : DEALLOCATE (ftau)
468 : END IF
469 :
470 6 : CALL auxbas_pw_pool%give_back_pw(rhoh_tot_gspace)
471 6 : CALL auxbas_pw_pool%give_back_pw(vhout_rspace)
472 6 : CALL auxbas_pw_pool%give_back_pw(vhout_gspace)
473 :
474 14 : DO ispin = 1, nspins
475 8 : CALL auxbas_pw_pool%give_back_pw(rhoh_r(ispin))
476 8 : CALL auxbas_pw_pool%give_back_pw(rhoh_g(ispin))
477 14 : CALL auxbas_pw_pool%give_back_pw(fhxc_rspace(ispin))
478 : END DO
479 6 : DEALLOCATE (rhoh_r, rhoh_g, fhxc_rspace)
480 :
481 6 : CALL timestop(handle)
482 :
483 18 : END SUBROUTINE harris_forces
484 :
485 : ! **************************************************************************************************
486 : !> \brief ...
487 : !> \param qs_env ...
488 : !> \param rhoin ...
489 : !> \param v_rspace ...
490 : !> \param calculate_forces ...
491 : ! **************************************************************************************************
492 6 : SUBROUTINE calculate_harris_integrals(qs_env, rhoin, v_rspace, calculate_forces)
493 : TYPE(qs_environment_type), POINTER :: qs_env
494 : TYPE(harris_rhoin_type), INTENT(INOUT) :: rhoin
495 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
496 : LOGICAL, INTENT(IN) :: calculate_forces
497 :
498 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_harris_integrals'
499 :
500 : INTEGER :: handle, i1, i2, iatom, ikind, ilocal, &
501 : ispin, n, nkind, nlocal, nspin
502 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: integral, vector
503 : TYPE(distribution_1d_type), POINTER :: local_particles
504 : TYPE(mp_para_env_type), POINTER :: para_env
505 :
506 6 : CALL timeset(routineN, handle)
507 :
508 6 : CALL get_qs_env(qs_env, para_env=para_env, local_particles=local_particles)
509 :
510 18 : ALLOCATE (vector(rhoin%nbas))
511 12 : ALLOCATE (integral(rhoin%nbas))
512 :
513 6 : nkind = SIZE(rhoin%rhovec, 1)
514 6 : nspin = SIZE(rhoin%rhovec, 2)
515 :
516 14 : DO ispin = 1, nspin
517 8 : vector = 0.0_dp
518 8 : integral = 0.0_dp
519 32 : DO ikind = 1, nkind
520 24 : nlocal = local_particles%n_el(ikind)
521 48 : DO ilocal = 1, nlocal
522 16 : iatom = local_particles%list(ikind)%array(ilocal)
523 16 : i1 = rhoin%basptr(iatom, 1)
524 16 : i2 = rhoin%basptr(iatom, 2)
525 16 : n = i2 - i1 + 1
526 144 : vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal)
527 : END DO
528 : END DO
529 8 : CALL para_env%sum(vector)
530 : !
531 : CALL integrate_function(qs_env, v_rspace(ispin), vector, integral, &
532 8 : calculate_forces, rhoin%basis_type)
533 38 : DO ikind = 1, nkind
534 24 : nlocal = local_particles%n_el(ikind)
535 48 : DO ilocal = 1, nlocal
536 16 : iatom = local_particles%list(ikind)%array(ilocal)
537 16 : i1 = rhoin%basptr(iatom, 1)
538 16 : i2 = rhoin%basptr(iatom, 2)
539 16 : n = i2 - i1 + 1
540 144 : rhoin%intvec(ikind, ispin)%rvecs(1:n, ilocal) = integral(i1:i2)
541 : END DO
542 : END DO
543 : END DO
544 :
545 6 : DEALLOCATE (vector, integral)
546 :
547 6 : CALL timestop(handle)
548 :
549 6 : END SUBROUTINE calculate_harris_integrals
550 :
551 : END MODULE qs_harris_methods
|
