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 Harris method environment setup and handling
10 : !> \par History
11 : !> 2024.07 created
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_harris_utils
15 : USE atom_kind_orbitals, ONLY: calculate_atomic_density
16 : USE atomic_kind_types, ONLY: atomic_kind_type
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE cell_types, ONLY: cell_type
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
22 : dbcsr_create,&
23 : dbcsr_p_type,&
24 : dbcsr_release,&
25 : dbcsr_set
26 : USE cp_log_handling, ONLY: cp_get_default_logger,&
27 : cp_logger_get_default_io_unit,&
28 : cp_logger_get_default_unit_nr,&
29 : cp_logger_type
30 : USE distribution_1d_types, ONLY: distribution_1d_type
31 : USE ec_methods, ONLY: create_kernel
32 : USE input_constants, ONLY: hden_atomic,&
33 : hfun_harris,&
34 : horb_default
35 : USE input_section_types, ONLY: section_vals_get_subs_vals,&
36 : section_vals_type,&
37 : section_vals_val_get
38 : USE kinds, ONLY: dp
39 : USE message_passing, ONLY: mp_para_env_type
40 : USE particle_types, ONLY: particle_type
41 : USE pw_env_types, ONLY: pw_env_get,&
42 : pw_env_type
43 : USE pw_grid_types, ONLY: pw_grid_type
44 : USE pw_methods, ONLY: pw_axpy,&
45 : pw_copy,&
46 : pw_integral_ab,&
47 : pw_integrate_function,&
48 : pw_scale,&
49 : pw_transfer,&
50 : pw_zero
51 : USE pw_poisson_methods, ONLY: pw_poisson_solve
52 : USE pw_poisson_types, ONLY: pw_poisson_type
53 : USE pw_pool_types, ONLY: pw_pool_type
54 : USE pw_types, ONLY: pw_c1d_gs_type,&
55 : pw_r3d_rs_type
56 : USE qs_collocate_density, ONLY: calculate_rho_elec,&
57 : collocate_function
58 : USE qs_energy_types, ONLY: qs_energy_type
59 : USE qs_environment_types, ONLY: get_qs_env,&
60 : qs_environment_type
61 : USE qs_force_types, ONLY: qs_force_type
62 : USE qs_harris_types, ONLY: harris_energy_type,&
63 : harris_print_energy,&
64 : harris_rhoin_type,&
65 : harris_type
66 : USE qs_integrate_potential, ONLY: integrate_function,&
67 : integrate_v_core_rspace,&
68 : integrate_v_rspace
69 : USE qs_kind_types, ONLY: get_qs_kind,&
70 : qs_kind_type
71 : USE qs_ks_types, ONLY: qs_ks_env_type
72 : USE qs_rho_types, ONLY: qs_rho_get,&
73 : qs_rho_type
74 : #include "./base/base_uses.f90"
75 :
76 : IMPLICIT NONE
77 :
78 : PRIVATE
79 :
80 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_harris_utils'
81 :
82 : PUBLIC :: harris_env_create, harris_write_input, harris_density_update, calculate_harris_density, &
83 : harris_energy_correction, harris_set_potentials
84 :
85 : CONTAINS
86 :
87 : ! **************************************************************************************************
88 : !> \brief Allocates and intitializes harris_env
89 : !> \param qs_env The QS environment
90 : !> \param harris_env The Harris method environment (the object to create)
91 : !> \param harris_section The Harris method input section
92 : !> \par History
93 : !> 2024.07 created
94 : !> \author JGH
95 : ! **************************************************************************************************
96 7404 : SUBROUTINE harris_env_create(qs_env, harris_env, harris_section)
97 : TYPE(qs_environment_type), POINTER :: qs_env
98 : TYPE(harris_type), POINTER :: harris_env
99 : TYPE(section_vals_type), OPTIONAL, POINTER :: harris_section
100 :
101 7404 : CPASSERT(.NOT. ASSOCIATED(harris_env))
102 7404 : ALLOCATE (harris_env)
103 7404 : CALL init_harris_env(qs_env, harris_env, harris_section)
104 :
105 7404 : END SUBROUTINE harris_env_create
106 :
107 : ! **************************************************************************************************
108 : !> \brief Initializes Harris method environment
109 : !> \param qs_env The QS environment
110 : !> \param harris_env The Harris method environment
111 : !> \param harris_section The Harris method input section
112 : !> \par History
113 : !> 2024.07 created
114 : !> \author JGH
115 : ! **************************************************************************************************
116 7404 : SUBROUTINE init_harris_env(qs_env, harris_env, harris_section)
117 : TYPE(qs_environment_type), POINTER :: qs_env
118 : TYPE(harris_type), POINTER :: harris_env
119 : TYPE(section_vals_type), OPTIONAL, POINTER :: harris_section
120 :
121 : CHARACTER(LEN=*), PARAMETER :: routineN = 'init_harris_env'
122 :
123 : INTEGER :: handle, unit_nr
124 : TYPE(cp_logger_type), POINTER :: logger
125 :
126 7404 : CALL timeset(routineN, handle)
127 :
128 7404 : IF (qs_env%harris_method) THEN
129 :
130 6 : CPASSERT(PRESENT(harris_section))
131 : ! get a useful output_unit
132 6 : logger => cp_get_default_logger()
133 6 : IF (logger%para_env%is_source()) THEN
134 3 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
135 : ELSE
136 : unit_nr = -1
137 : END IF
138 :
139 : CALL section_vals_val_get(harris_section, "ENERGY_FUNCTIONAL", &
140 6 : i_val=harris_env%energy_functional)
141 : CALL section_vals_val_get(harris_section, "DENSITY_SOURCE", &
142 6 : i_val=harris_env%density_source)
143 : CALL section_vals_val_get(harris_section, "ORBITAL_BASIS", &
144 6 : i_val=harris_env%orbital_basis)
145 : !
146 : CALL section_vals_val_get(harris_section, "DEBUG_FORCES", &
147 6 : l_val=harris_env%debug_forces)
148 : CALL section_vals_val_get(harris_section, "DEBUG_STRESS", &
149 6 : l_val=harris_env%debug_stress)
150 :
151 : END IF
152 :
153 7404 : CALL timestop(handle)
154 :
155 7404 : END SUBROUTINE init_harris_env
156 :
157 : ! **************************************************************************************************
158 : !> \brief Print out the Harris method input section
159 : !>
160 : !> \param harris_env ...
161 : !> \par History
162 : !> 2024.07 created [JGH]
163 : !> \author JGH
164 : ! **************************************************************************************************
165 6 : SUBROUTINE harris_write_input(harris_env)
166 : TYPE(harris_type), POINTER :: harris_env
167 :
168 : CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_write_input'
169 :
170 : INTEGER :: handle, unit_nr
171 : TYPE(cp_logger_type), POINTER :: logger
172 :
173 6 : CALL timeset(routineN, handle)
174 :
175 6 : logger => cp_get_default_logger()
176 6 : IF (logger%para_env%is_source()) THEN
177 3 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
178 : ELSE
179 : unit_nr = -1
180 : END IF
181 :
182 3 : IF (unit_nr > 0) THEN
183 :
184 : WRITE (unit_nr, '(/,T2,A)') &
185 3 : "!"//REPEAT("-", 29)//" Harris Model "//REPEAT("-", 29)//"!"
186 :
187 : ! Type of energy functional
188 6 : SELECT CASE (harris_env%energy_functional)
189 : CASE (hfun_harris)
190 3 : WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Functional: ", "Harris"
191 : END SELECT
192 : ! density source
193 6 : SELECT CASE (harris_env%density_source)
194 : CASE (hden_atomic)
195 3 : WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model density: Type", " Atomic kind density"
196 : END SELECT
197 3 : WRITE (unit_nr, '(T2,A,T71,A10)') "Harris model density: Basis type", &
198 6 : ADJUSTR(TRIM(harris_env%rhoin%basis_type))
199 3 : WRITE (unit_nr, '(T2,A,T71,I10)') "Harris model density: Number of basis functions", &
200 6 : harris_env%rhoin%nbas
201 : ! orbital basis
202 6 : SELECT CASE (harris_env%orbital_basis)
203 : CASE (horb_default)
204 3 : WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model basis: ", "Atomic kind orbitals"
205 : END SELECT
206 :
207 3 : WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
208 3 : WRITE (unit_nr, '()')
209 :
210 : END IF ! unit_nr
211 :
212 6 : CALL timestop(handle)
213 :
214 6 : END SUBROUTINE harris_write_input
215 :
216 : ! **************************************************************************************************
217 : !> \brief ...
218 : !> \param qs_env ...
219 : !> \param harris_env ...
220 : ! **************************************************************************************************
221 30 : SUBROUTINE harris_density_update(qs_env, harris_env)
222 : TYPE(qs_environment_type), POINTER :: qs_env
223 : TYPE(harris_type), POINTER :: harris_env
224 :
225 : CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_density_update'
226 :
227 : INTEGER :: handle, i, ikind, ngto, nkind, nset, nsgf
228 30 : INTEGER, DIMENSION(:), POINTER :: lmax, npgf
229 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coef
230 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: density
231 30 : REAL(KIND=dp), DIMENSION(:), POINTER :: norm
232 30 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
233 30 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
234 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
235 : TYPE(atomic_kind_type), POINTER :: atomic_kind
236 : TYPE(gto_basis_set_type), POINTER :: basis_set
237 30 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
238 : TYPE(qs_kind_type), POINTER :: qs_kind
239 :
240 30 : CALL timeset(routineN, handle)
241 :
242 60 : SELECT CASE (harris_env%density_source)
243 : CASE (hden_atomic)
244 30 : IF (.NOT. harris_env%rhoin%frozen) THEN
245 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
246 6 : nkind=nkind)
247 22 : DO ikind = 1, nkind
248 16 : atomic_kind => atomic_kind_set(ikind)
249 16 : qs_kind => qs_kind_set(ikind)
250 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, &
251 16 : basis_type=harris_env%rhoin%basis_type)
252 : CALL get_gto_basis_set(gto_basis_set=basis_set, nset=nset, lmax=lmax, nsgf=nsgf, &
253 16 : npgf=npgf, norm_cgf=norm, zet=zet, gcc=gcc)
254 16 : IF (nset /= 1 .OR. lmax(1) /= 0 .OR. npgf(1) /= nsgf) THEN
255 0 : CPABORT("RHOIN illegal basis type")
256 : END IF
257 122 : DO i = 1, npgf(1)
258 1746 : IF (SUM(ABS(gcc(1:npgf(1), i, 1))) /= MAXVAL(ABS(gcc(1:npgf(1), i, 1)))) THEN
259 0 : CPABORT("RHOIN illegal basis type")
260 : END IF
261 : END DO
262 : !
263 16 : ngto = npgf(1)
264 48 : ALLOCATE (density(ngto, 2))
265 122 : density(1:ngto, 1) = zet(1:ngto, 1)
266 122 : density(1:ngto, 2) = 0.0_dp
267 : CALL calculate_atomic_density(density, atomic_kind, qs_kind, ngto, &
268 16 : optbasis=.FALSE., confine=.TRUE.)
269 48 : ALLOCATE (coef(ngto))
270 122 : DO i = 1, ngto
271 122 : coef(i) = density(i, 2)/gcc(i, i, 1)/norm(i)
272 : END DO
273 16 : IF (harris_env%rhoin%nspin == 2) THEN
274 0 : DO i = 1, SIZE(harris_env%rhoin%rhovec(ikind, 1)%rvecs, 2)
275 0 : harris_env%rhoin%rhovec(ikind, 1)%rvecs(1:ngto, i) = coef(1:ngto)*0.5_dp
276 0 : harris_env%rhoin%rhovec(ikind, 2)%rvecs(1:ngto, i) = coef(1:ngto)*0.5_dp
277 : END DO
278 : ELSE
279 30 : DO i = 1, SIZE(harris_env%rhoin%rhovec(ikind, 1)%rvecs, 2)
280 120 : harris_env%rhoin%rhovec(ikind, 1)%rvecs(1:ngto, i) = coef(1:ngto)
281 : END DO
282 : END IF
283 38 : DEALLOCATE (density, coef)
284 : END DO
285 6 : harris_env%rhoin%frozen = .TRUE.
286 : END IF
287 : CASE DEFAULT
288 30 : CPABORT("Illeagal value of harris_env%density_source")
289 : END SELECT
290 :
291 30 : CALL timestop(handle)
292 :
293 60 : END SUBROUTINE harris_density_update
294 :
295 : ! **************************************************************************************************
296 : !> \brief ...
297 : !> \param qs_env ...
298 : !> \param rhoin ...
299 : !> \param rho_struct ...
300 : ! **************************************************************************************************
301 42 : SUBROUTINE calculate_harris_density(qs_env, rhoin, rho_struct)
302 : TYPE(qs_environment_type), POINTER :: qs_env
303 : TYPE(harris_rhoin_type), INTENT(IN) :: rhoin
304 : TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
305 :
306 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_harris_density'
307 :
308 : INTEGER :: handle, i1, i2, iatom, ikind, ilocal, &
309 : ispin, n, nkind, nlocal, nspin
310 : REAL(KIND=dp) :: eps_rho_rspace
311 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vector
312 42 : REAL(KIND=dp), DIMENSION(:), POINTER :: total_rho
313 42 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
314 : TYPE(cell_type), POINTER :: cell
315 : TYPE(dft_control_type), POINTER :: dft_control
316 : TYPE(distribution_1d_type), POINTER :: local_particles
317 : TYPE(mp_para_env_type), POINTER :: para_env
318 42 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
319 42 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_gspace
320 : TYPE(pw_env_type), POINTER :: pw_env
321 42 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_rspace
322 42 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
323 :
324 42 : CALL timeset(routineN, handle)
325 :
326 42 : CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
327 42 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
328 : CALL get_qs_env(qs_env, &
329 : atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
330 : local_particles=local_particles, &
331 42 : qs_kind_set=qs_kind_set, cell=cell, pw_env=pw_env)
332 :
333 : CALL qs_rho_get(rho_struct, rho_r=rho_rspace, rho_g=rho_gspace, &
334 42 : tot_rho_r=total_rho)
335 :
336 126 : ALLOCATE (vector(rhoin%nbas))
337 :
338 42 : nkind = SIZE(rhoin%rhovec, 1)
339 42 : nspin = SIZE(rhoin%rhovec, 2)
340 :
341 84 : DO ispin = 1, nspin
342 1158 : vector = 0.0_dp
343 166 : DO ikind = 1, nkind
344 124 : nlocal = local_particles%n_el(ikind)
345 252 : DO ilocal = 1, nlocal
346 86 : iatom = local_particles%list(ikind)%array(ilocal)
347 86 : i1 = rhoin%basptr(iatom, 1)
348 86 : i2 = rhoin%basptr(iatom, 2)
349 86 : n = i2 - i1 + 1
350 768 : vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal)
351 : END DO
352 : END DO
353 42 : CALL para_env%sum(vector)
354 : !
355 : CALL collocate_function(vector, rho_rspace(ispin), rho_gspace(ispin), &
356 : atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
357 42 : eps_rho_rspace, rhoin%basis_type)
358 84 : total_rho(ispin) = pw_integrate_function(rho_rspace(ispin), isign=-1)
359 : END DO
360 :
361 42 : DEALLOCATE (vector)
362 :
363 42 : CALL timestop(handle)
364 :
365 42 : END SUBROUTINE calculate_harris_density
366 :
367 : ! **************************************************************************************************
368 : !> \brief ...
369 : !> \param qs_env ...
370 : !> \param rhoin ...
371 : !> \param v_rspace ...
372 : !> \param calculate_forces ...
373 : ! **************************************************************************************************
374 4 : SUBROUTINE calculate_harris_integrals(qs_env, rhoin, v_rspace, calculate_forces)
375 : TYPE(qs_environment_type), POINTER :: qs_env
376 : TYPE(harris_rhoin_type), INTENT(INOUT) :: rhoin
377 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
378 : LOGICAL, INTENT(IN) :: calculate_forces
379 :
380 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_harris_integrals'
381 :
382 : INTEGER :: handle, i1, i2, iatom, ikind, ilocal, &
383 : ispin, n, nkind, nlocal, nspin
384 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: integral, vector
385 : TYPE(distribution_1d_type), POINTER :: local_particles
386 : TYPE(mp_para_env_type), POINTER :: para_env
387 :
388 4 : CALL timeset(routineN, handle)
389 :
390 4 : CALL get_qs_env(qs_env, para_env=para_env, local_particles=local_particles)
391 :
392 12 : ALLOCATE (vector(rhoin%nbas))
393 12 : ALLOCATE (integral(rhoin%nbas))
394 :
395 4 : nkind = SIZE(rhoin%rhovec, 1)
396 4 : nspin = SIZE(rhoin%rhovec, 2)
397 :
398 8 : DO ispin = 1, nspin
399 108 : vector = 0.0_dp
400 108 : integral = 0.0_dp
401 16 : DO ikind = 1, nkind
402 12 : nlocal = local_particles%n_el(ikind)
403 24 : DO ilocal = 1, nlocal
404 8 : iatom = local_particles%list(ikind)%array(ilocal)
405 8 : i1 = rhoin%basptr(iatom, 1)
406 8 : i2 = rhoin%basptr(iatom, 2)
407 8 : n = i2 - i1 + 1
408 72 : vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal)
409 : END DO
410 : END DO
411 4 : CALL para_env%sum(vector)
412 : !
413 : CALL integrate_function(qs_env, v_rspace(ispin), vector, integral, &
414 4 : calculate_forces, rhoin%basis_type)
415 20 : DO ikind = 1, nkind
416 12 : nlocal = local_particles%n_el(ikind)
417 24 : DO ilocal = 1, nlocal
418 8 : iatom = local_particles%list(ikind)%array(ilocal)
419 8 : i1 = rhoin%basptr(iatom, 1)
420 8 : i2 = rhoin%basptr(iatom, 2)
421 8 : n = i2 - i1 + 1
422 72 : rhoin%intvec(ikind, ispin)%rvecs(1:n, ilocal) = integral(i1:i2)
423 : END DO
424 : END DO
425 : END DO
426 :
427 4 : DEALLOCATE (vector, integral)
428 :
429 4 : CALL timestop(handle)
430 :
431 4 : END SUBROUTINE calculate_harris_integrals
432 :
433 : ! **************************************************************************************************
434 : !> \brief ...
435 : !> \param harris_env ...
436 : !> \param vh_rspace ...
437 : !> \param vxc_rspace ...
438 : ! **************************************************************************************************
439 54 : SUBROUTINE harris_set_potentials(harris_env, vh_rspace, vxc_rspace)
440 : TYPE(harris_type), POINTER :: harris_env
441 : TYPE(pw_r3d_rs_type), INTENT(IN) :: vh_rspace
442 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace
443 :
444 : INTEGER :: iab, ispin, nspins
445 : TYPE(pw_grid_type), POINTER :: pw_grid
446 :
447 : ! release possible old potentials
448 54 : IF (ASSOCIATED(harris_env%vh_rspace%pw_grid)) THEN
449 48 : CALL harris_env%vh_rspace%release()
450 : END IF
451 54 : IF (ASSOCIATED(harris_env%vxc_rspace)) THEN
452 96 : DO iab = 1, SIZE(harris_env%vxc_rspace)
453 96 : CALL harris_env%vxc_rspace(iab)%release()
454 : END DO
455 48 : DEALLOCATE (harris_env%vxc_rspace)
456 : END IF
457 :
458 : ! generate new potential data structures
459 54 : nspins = harris_env%rhoin%nspin
460 216 : ALLOCATE (harris_env%vxc_rspace(nspins))
461 :
462 54 : pw_grid => vh_rspace%pw_grid
463 54 : CALL harris_env%vh_rspace%create(pw_grid)
464 108 : DO ispin = 1, nspins
465 108 : CALL harris_env%vxc_rspace(ispin)%create(pw_grid)
466 : END DO
467 :
468 : ! copy potentials
469 54 : CALL pw_transfer(vh_rspace, harris_env%vh_rspace)
470 54 : IF (ASSOCIATED(vxc_rspace)) THEN
471 108 : DO ispin = 1, nspins
472 54 : CALL pw_transfer(vxc_rspace(ispin), harris_env%vxc_rspace(ispin))
473 108 : CALL pw_scale(harris_env%vxc_rspace(ispin), vxc_rspace(ispin)%pw_grid%dvol)
474 : END DO
475 : ELSE
476 0 : DO ispin = 1, nspins
477 0 : CALL pw_zero(harris_env%vxc_rspace(ispin))
478 : END DO
479 : END IF
480 :
481 54 : END SUBROUTINE harris_set_potentials
482 :
483 : ! **************************************************************************************************
484 : !> \brief ...
485 : !> \param qs_env ...
486 : !> \param calculate_forces ...
487 : ! **************************************************************************************************
488 30 : SUBROUTINE harris_energy_correction(qs_env, calculate_forces)
489 : TYPE(qs_environment_type), POINTER :: qs_env
490 : LOGICAL, INTENT(IN) :: calculate_forces
491 :
492 : CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_energy_correction'
493 :
494 : INTEGER :: handle, iounit, ispin, nspins
495 : REAL(KIND=dp) :: dvol, ec, eh, exc, vxc
496 : TYPE(cp_logger_type), POINTER :: logger
497 : TYPE(harris_energy_type), POINTER :: energy
498 : TYPE(harris_type), POINTER :: harris_env
499 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
500 : TYPE(pw_env_type), POINTER :: pw_env
501 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
502 : TYPE(pw_r3d_rs_type) :: core_rspace
503 30 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
504 : TYPE(qs_energy_type), POINTER :: ks_energy
505 : TYPE(qs_rho_type), POINTER :: rho
506 :
507 : MARK_USED(calculate_forces)
508 :
509 30 : CALL timeset(routineN, handle)
510 :
511 30 : CALL get_qs_env(qs_env, harris_env=harris_env, energy=ks_energy)
512 30 : energy => harris_env%energy
513 30 : energy%eband = ks_energy%band
514 30 : energy%ewald_correction = ks_energy%core_overlap + ks_energy%core_self
515 30 : energy%dispersion = ks_energy%dispersion
516 :
517 30 : nspins = harris_env%rhoin%nspin
518 :
519 30 : CALL get_qs_env(qs_env, rho=rho, rho_core=rho_core)
520 30 : CALL qs_rho_get(rho, rho_r=rho_r)
521 :
522 30 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
523 30 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
524 30 : CALL auxbas_pw_pool%create_pw(core_rspace)
525 30 : CALL pw_transfer(rho_core, core_rspace)
526 :
527 30 : dvol = harris_env%vh_rspace%pw_grid%dvol
528 30 : eh = 0.0_dp
529 60 : DO ispin = 1, nspins
530 60 : eh = eh + pw_integral_ab(rho_r(ispin), harris_env%vh_rspace)/dvol
531 : END DO
532 30 : ec = pw_integral_ab(core_rspace, harris_env%vh_rspace)/dvol
533 30 : eh = 0.5_dp*(eh + ec)
534 30 : energy%eh_correction = ec - eh
535 :
536 30 : exc = ks_energy%exc
537 30 : vxc = 0.0_dp
538 30 : IF (ASSOCIATED(harris_env%vxc_rspace)) THEN
539 60 : DO ispin = 1, nspins
540 : vxc = vxc + pw_integral_ab(rho_r(ispin), harris_env%vxc_rspace(ispin))/ &
541 60 : harris_env%vxc_rspace(ispin)%pw_grid%dvol
542 : END DO
543 : END IF
544 30 : energy%exc_correction = exc - vxc
545 :
546 : ! Total Harris model energy
547 : energy%eharris = energy%eband + energy%eh_correction + energy%exc_correction + &
548 30 : energy%ewald_correction + energy%dispersion
549 :
550 30 : CALL auxbas_pw_pool%give_back_pw(core_rspace)
551 :
552 30 : ks_energy%total = ks_energy%total + ks_energy%core
553 30 : ks_energy%nonscf_correction = energy%eharris - ks_energy%total
554 30 : ks_energy%total = energy%eharris
555 :
556 30 : logger => cp_get_default_logger()
557 30 : iounit = cp_logger_get_default_io_unit(logger)
558 :
559 30 : CALL harris_print_energy(iounit, energy)
560 :
561 30 : IF (calculate_forces) THEN
562 4 : CALL harris_forces(qs_env, iounit)
563 : END IF
564 :
565 30 : CALL timestop(handle)
566 :
567 30 : END SUBROUTINE harris_energy_correction
568 :
569 : ! **************************************************************************************************
570 : !> \brief ...
571 : !> \param qs_env ...
572 : !> \param iounit ...
573 : ! **************************************************************************************************
574 4 : SUBROUTINE harris_forces(qs_env, iounit)
575 : TYPE(qs_environment_type), POINTER :: qs_env
576 : INTEGER, INTENT(IN) :: iounit
577 :
578 : CHARACTER(LEN=*), PARAMETER :: routineN = 'harris_forces'
579 : LOGICAL, PARAMETER :: debug_forces = .TRUE.
580 :
581 : INTEGER :: handle, ispin, nspins
582 : REAL(KIND=dp) :: ehartree
583 : REAL(KIND=dp), DIMENSION(3) :: fodeb
584 : TYPE(dbcsr_p_type) :: scrm
585 4 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: rhoh_ao, smat
586 : TYPE(harris_type), POINTER :: harris_env
587 : TYPE(mp_para_env_type), POINTER :: para_env
588 : TYPE(pw_c1d_gs_type) :: rhoh_tot_gspace, vhout_gspace
589 4 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g, rhoh_g
590 : TYPE(pw_c1d_gs_type), POINTER :: rho_core
591 : TYPE(pw_env_type), POINTER :: pw_env
592 : TYPE(pw_poisson_type), POINTER :: poisson_env
593 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
594 : TYPE(pw_r3d_rs_type) :: vhout_rspace, vhxc_rspace
595 4 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: fhxc_rspace, ftau, fxc, rho_r, rhoh_r, &
596 4 : tauh_r
597 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
598 : TYPE(qs_ks_env_type), POINTER :: ks_env
599 : TYPE(qs_rho_type), POINTER :: rho
600 : TYPE(section_vals_type), POINTER :: xc_section
601 :
602 4 : CALL timeset(routineN, handle)
603 :
604 : IF (debug_forces) THEN
605 4 : IF (iounit > 0) WRITE (iounit, "(/,T3,A)") &
606 2 : "DEBUG:: Harris Method Forces (density dependent)"
607 : END IF
608 :
609 4 : CALL get_qs_env(qs_env, harris_env=harris_env, force=force, para_env=para_env)
610 4 : nspins = harris_env%rhoin%nspin
611 :
612 4 : CALL get_qs_env(qs_env, rho=rho, rho_core=rho_core, matrix_s=smat)
613 : ! Warning: rho_ao = output DM; rho_r = rhoin
614 4 : CALL qs_rho_get(rho, rho_ao=rhoh_ao, rho_r=rho_r, rho_g=rho_g)
615 4 : ALLOCATE (scrm%matrix)
616 4 : CALL dbcsr_create(scrm%matrix, template=rhoh_ao(1)%matrix)
617 4 : CALL dbcsr_copy(scrm%matrix, smat(1)%matrix)
618 4 : CALL dbcsr_set(scrm%matrix, 0.0_dp)
619 :
620 4 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, ks_env=ks_env)
621 4 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
622 4 : CALL auxbas_pw_pool%create_pw(vhxc_rspace)
623 :
624 16 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
625 8 : DO ispin = 1, nspins
626 4 : CALL pw_copy(harris_env%vh_rspace, vhxc_rspace)
627 4 : CALL pw_axpy(harris_env%vxc_rspace(ispin), vhxc_rspace)
628 : CALL integrate_v_rspace(v_rspace=vhxc_rspace, &
629 : hmat=scrm, pmat=rhoh_ao(ispin), &
630 8 : qs_env=qs_env, calculate_forces=.TRUE.)
631 : END DO
632 : IF (debug_forces) THEN
633 16 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
634 4 : CALL para_env%sum(fodeb)
635 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: P*(Vh[in]+Vxc)", fodeb
636 : END IF
637 :
638 4 : CALL dbcsr_release(scrm%matrix)
639 4 : DEALLOCATE (scrm%matrix)
640 4 : CALL auxbas_pw_pool%give_back_pw(vhxc_rspace)
641 :
642 28 : ALLOCATE (rhoh_r(nspins), rhoh_g(nspins))
643 8 : DO ispin = 1, nspins
644 4 : CALL auxbas_pw_pool%create_pw(rhoh_r(ispin))
645 8 : CALL auxbas_pw_pool%create_pw(rhoh_g(ispin))
646 : END DO
647 4 : CALL auxbas_pw_pool%create_pw(rhoh_tot_gspace)
648 4 : CALL pw_copy(rho_core, rhoh_tot_gspace)
649 8 : DO ispin = 1, nspins
650 : CALL calculate_rho_elec(ks_env=ks_env, matrix_p=rhoh_ao(ispin)%matrix, &
651 4 : rho=rhoh_r(ispin), rho_gspace=rhoh_g(ispin))
652 8 : CALL pw_axpy(rhoh_g(ispin), rhoh_tot_gspace)
653 : END DO
654 : ! no meta functionals here
655 4 : NULLIFY (tauh_r)
656 :
657 4 : CALL auxbas_pw_pool%create_pw(vhout_rspace)
658 4 : CALL auxbas_pw_pool%create_pw(vhout_gspace)
659 4 : CALL pw_env_get(pw_env, poisson_env=poisson_env)
660 : !
661 4 : CALL pw_poisson_solve(poisson_env, rhoh_tot_gspace, ehartree, vhout_gspace)
662 : !
663 4 : CALL pw_transfer(vhout_gspace, vhout_rspace)
664 4 : CALL pw_scale(vhout_rspace, vhout_rspace%pw_grid%dvol)
665 :
666 16 : IF (debug_forces) fodeb(1:3) = force(1)%rho_core(1:3, 1)
667 4 : CALL integrate_v_core_rspace(vhout_rspace, qs_env)
668 : IF (debug_forces) THEN
669 16 : fodeb(1:3) = force(1)%rho_core(1:3, 1) - fodeb(1:3)
670 4 : CALL para_env%sum(fodeb)
671 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: Vh[out]*dncore ", fodeb
672 : END IF
673 :
674 12 : ALLOCATE (fhxc_rspace(nspins))
675 8 : DO ispin = 1, nspins
676 8 : CALL auxbas_pw_pool%create_pw(fhxc_rspace(ispin))
677 : END DO
678 : ! vh = vh[out] - vh[in]
679 4 : CALL pw_axpy(harris_env%vh_rspace, vhout_rspace, alpha=-1._dp, beta=1.0_dp)
680 : ! kernel fxc
681 : ! drho = rho[out] - rho[in]
682 8 : DO ispin = 1, nspins
683 4 : CALL pw_axpy(rho_r(ispin), rhoh_r(ispin), alpha=-1._dp, beta=1.0_dp)
684 8 : CALL pw_axpy(rho_g(ispin), rhoh_g(ispin), alpha=-1._dp, beta=1.0_dp)
685 : END DO
686 4 : xc_section => section_vals_get_subs_vals(qs_env%input, "DFT%XC")
687 4 : NULLIFY (fxc, ftau)
688 : CALL create_kernel(qs_env, vxc=fxc, vxc_tau=ftau, &
689 : rho=rho, rho1_r=rhoh_r, rho1_g=rhoh_g, tau1_r=tauh_r, &
690 4 : xc_section=xc_section)
691 4 : CPASSERT(.NOT. ASSOCIATED(ftau))
692 :
693 8 : DO ispin = 1, nspins
694 4 : CALL pw_copy(vhout_rspace, fhxc_rspace(ispin))
695 8 : IF (ASSOCIATED(fxc)) THEN
696 4 : CALL pw_scale(fxc(ispin), fxc(ispin)%pw_grid%dvol)
697 4 : CALL pw_axpy(fxc(ispin), fhxc_rspace(ispin))
698 : END IF
699 : END DO
700 :
701 16 : IF (debug_forces) fodeb(1:3) = force(1)%rho_elec(1:3, 1)
702 4 : CALL calculate_harris_integrals(qs_env, harris_env%rhoin, fhxc_rspace, .TRUE.)
703 : IF (debug_forces) THEN
704 16 : fodeb(1:3) = force(1)%rho_elec(1:3, 1) - fodeb(1:3)
705 4 : CALL para_env%sum(fodeb)
706 4 : IF (iounit > 0) WRITE (iounit, "(T3,A,T33,3F16.8)") "DEBUG:: (dVh+fxc)*dn[in] ", fodeb
707 : END IF
708 :
709 4 : IF (ASSOCIATED(fxc)) THEN
710 8 : DO ispin = 1, nspins
711 8 : CALL auxbas_pw_pool%give_back_pw(fxc(ispin))
712 : END DO
713 4 : DEALLOCATE (fxc)
714 : END IF
715 4 : IF (ASSOCIATED(ftau)) THEN
716 0 : DO ispin = 1, nspins
717 0 : CALL auxbas_pw_pool%give_back_pw(ftau(ispin))
718 : END DO
719 0 : DEALLOCATE (ftau)
720 : END IF
721 :
722 4 : CALL auxbas_pw_pool%give_back_pw(rhoh_tot_gspace)
723 4 : CALL auxbas_pw_pool%give_back_pw(vhout_rspace)
724 4 : CALL auxbas_pw_pool%give_back_pw(vhout_gspace)
725 :
726 8 : DO ispin = 1, nspins
727 4 : CALL auxbas_pw_pool%give_back_pw(rhoh_r(ispin))
728 4 : CALL auxbas_pw_pool%give_back_pw(rhoh_g(ispin))
729 8 : CALL auxbas_pw_pool%give_back_pw(fhxc_rspace(ispin))
730 : END DO
731 4 : DEALLOCATE (rhoh_r, rhoh_g, fhxc_rspace)
732 :
733 4 : CALL timestop(handle)
734 :
735 8 : END SUBROUTINE harris_forces
736 :
737 : END MODULE qs_harris_utils
|
