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