(git:591cf04)
Loading...
Searching...
No Matches
qs_harris_utils.F
Go to the documentation of this file.
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 environment setup and handling
10!> \par History
11!> 2024.07 created
12!> \author JGH
13! **************************************************************************************************
14MODULE 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_log_handling, ONLY: cp_get_default_logger,&
22 cp_logger_get_default_unit_nr,&
23 cp_logger_type
24 USE distribution_1d_types, ONLY: distribution_1d_type
25 USE input_constants, ONLY: hden_atomic,&
26 hden_cube,&
27 hden_cube_fit,&
28 hfit_least_squares,&
29 hfit_relative_entropy,&
30 hfun_harris,&
31 horb_default
32 USE input_section_types, ONLY: section_vals_type,&
33 section_vals_val_get
34 USE kinds, ONLY: dp
35 USE message_passing, ONLY: mp_para_env_type
36 USE particle_types, ONLY: particle_type
37 USE pw_env_types, ONLY: pw_env_type
38 USE pw_grid_types, ONLY: pw_grid_type
39 USE pw_methods, ONLY: pw_copy,&
40 pw_integrate_function,&
41 pw_scale,&
42 pw_transfer,&
43 pw_zero
44 USE pw_types, ONLY: pw_c1d_gs_type,&
45 pw_r3d_rs_type
46 USE qs_collocate_density, ONLY: collocate_function
47 USE qs_density_fit, ONLY: fit_constrained_density
48 USE qs_environment_types, ONLY: get_qs_env,&
49 qs_environment_type
50 USE qs_external_density, ONLY: read_cube_density
51 USE qs_harris_types, ONLY: harris_rhoin_type,&
52 harris_type
53 USE qs_integrate_potential, ONLY: integrate_function
54 USE qs_kind_types, ONLY: get_qs_kind,&
55 qs_kind_type
56 USE qs_rho_types, ONLY: qs_rho_get,&
57 qs_rho_type
58#include "./base/base_uses.f90"
59
60 IMPLICIT NONE
61
62 PRIVATE
63
64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_harris_utils'
65
66 PUBLIC :: harris_env_create, harris_write_input, harris_density_update, calculate_harris_density, &
67 harris_set_potentials
68
69CONTAINS
70
71! **************************************************************************************************
72!> \brief Allocates and intitializes harris_env
73!> \param qs_env The QS environment
74!> \param harris_env The Harris method environment (the object to create)
75!> \param harris_section The Harris method input section
76!> \par History
77!> 2024.07 created
78!> \author JGH
79! **************************************************************************************************
80 SUBROUTINE harris_env_create(qs_env, harris_env, harris_section)
81 TYPE(qs_environment_type), POINTER :: qs_env
82 TYPE(harris_type), POINTER :: harris_env
83 TYPE(section_vals_type), OPTIONAL, POINTER :: harris_section
84
85 cpassert(.NOT. ASSOCIATED(harris_env))
86 ALLOCATE (harris_env)
87 CALL init_harris_env(qs_env, harris_env, harris_section)
88
89 END SUBROUTINE harris_env_create
90
91! **************************************************************************************************
92!> \brief Initializes Harris method environment
93!> \param qs_env The QS environment
94!> \param harris_env The Harris method environment
95!> \param harris_section The Harris method input section
96!> \par History
97!> 2024.07 created
98!> \author JGH
99! **************************************************************************************************
100 SUBROUTINE init_harris_env(qs_env, harris_env, harris_section)
101 TYPE(qs_environment_type), POINTER :: qs_env
102 TYPE(harris_type), POINTER :: harris_env
103 TYPE(section_vals_type), OPTIONAL, POINTER :: harris_section
104
105 CHARACTER(LEN=*), PARAMETER :: routinen = 'init_harris_env'
106
107 INTEGER :: handle, unit_nr
108 TYPE(cp_logger_type), POINTER :: logger
109
110 CALL timeset(routinen, handle)
111
112 IF (qs_env%harris_method) THEN
113
114 cpassert(PRESENT(harris_section))
115 ! get a useful output_unit
116 logger => cp_get_default_logger()
117 IF (logger%para_env%is_source()) THEN
118 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
119 ELSE
120 unit_nr = -1
121 END IF
122
123 CALL section_vals_val_get(harris_section, "ENERGY_FUNCTIONAL", &
124 i_val=harris_env%energy_functional)
125 CALL section_vals_val_get(harris_section, "DENSITY_SOURCE", &
126 i_val=harris_env%density_source)
127 CALL section_vals_val_get(harris_section, "FILE_DENSITY", &
128 c_val=harris_env%density_filename)
129 CALL section_vals_val_get(harris_section, "FIT_MAX_ITER", &
130 i_val=harris_env%fit_max_iter)
131 CALL section_vals_val_get(harris_section, "FIT_METHOD", &
132 i_val=harris_env%fit_method)
133 CALL section_vals_val_get(harris_section, "FIT_EPS", &
134 r_val=harris_env%fit_eps)
135 CALL section_vals_val_get(harris_section, "FIT_STEP_SIZE", &
136 r_val=harris_env%fit_step_size)
137 CALL section_vals_val_get(harris_section, "FIT_MAX_BACKTRACK", &
138 i_val=harris_env%fit_max_backtrack)
139 CALL section_vals_val_get(harris_section, "FIT_TEMPERATURE", &
140 r_val=harris_env%fit_temperature)
141 CALL section_vals_val_get(harris_section, "FIT_RELATIVE_ENTROPY_WEIGHT", &
142 r_val=harris_env%fit_relative_entropy_weight)
143 CALL section_vals_val_get(harris_section, "DIRECT_DENSITY_MATRIX_ENERGY", &
144 l_val=harris_env%direct_density_matrix_energy)
145 CALL section_vals_val_get(harris_section, "ORBITAL_BASIS", &
146 i_val=harris_env%orbital_basis)
147 !
148 CALL section_vals_val_get(harris_section, "DEBUG_FORCES", &
149 l_val=harris_env%debug_forces)
150 CALL section_vals_val_get(harris_section, "DEBUG_STRESS", &
151 l_val=harris_env%debug_stress)
152
153 END IF
154
155 CALL timestop(handle)
156
157 END SUBROUTINE init_harris_env
158
159! **************************************************************************************************
160!> \brief Print out the Harris method input section
161!>
162!> \param harris_env ...
163!> \par History
164!> 2024.07 created [JGH]
165!> \author JGH
166! **************************************************************************************************
167 SUBROUTINE harris_write_input(harris_env)
168 TYPE(harris_type), POINTER :: harris_env
169
170 CHARACTER(LEN=*), PARAMETER :: routinen = 'harris_write_input'
171
172 INTEGER :: handle, unit_nr
173 TYPE(cp_logger_type), POINTER :: logger
174
175 CALL timeset(routinen, handle)
176
177 logger => cp_get_default_logger()
178 IF (logger%para_env%is_source()) THEN
179 unit_nr = cp_logger_get_default_unit_nr(logger, local=.true.)
180 ELSE
181 unit_nr = -1
182 END IF
183
184 IF (unit_nr > 0) THEN
185
186 WRITE (unit_nr, '(/,T2,A)') &
187 "!"//repeat("-", 29)//" Harris Model "//repeat("-", 29)//"!"
188
189 ! Type of energy functional
190 SELECT CASE (harris_env%energy_functional)
191 CASE (hfun_harris)
192 WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Functional: ", "Harris"
193 END SELECT
194 ! density source
195 SELECT CASE (harris_env%density_source)
196 CASE (hden_atomic)
197 WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model density: Type", " Atomic kind density"
198 CASE (hden_cube)
199 WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model density: Type", "Cube file"
200 WRITE (unit_nr, '(T2,A,T31,A)') "Harris model density: File", &
201 trim(harris_env%density_filename)
202 CASE (hden_cube_fit)
203 WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model density: Type", "Constrained cube fit"
204 WRITE (unit_nr, '(T2,A,T31,A)') "Harris model density: File", &
205 trim(harris_env%density_filename)
206 WRITE (unit_nr, '(T2,A,T61,I20)') "Harris density fit: Maximum iterations", &
207 harris_env%fit_max_iter
208 WRITE (unit_nr, '(T2,A,T61,ES20.8)') "Harris density fit: RMS target", &
209 harris_env%fit_eps
210 WRITE (unit_nr, '(T2,A,T61,ES20.8)') "Harris density fit: Initial step size", &
211 harris_env%fit_step_size
212 SELECT CASE (harris_env%fit_method)
213 CASE (hfit_least_squares)
214 WRITE (unit_nr, '(T2,A,T61,A20)') "Harris density fit: Objective", "Least squares"
215 CASE (hfit_relative_entropy)
216 WRITE (unit_nr, '(T2,A,T61,A20)') "Harris density fit: Objective", "Relative entropy"
217 WRITE (unit_nr, '(T2,A,T61,ES20.8)') "Harris density fit: Temperature", &
218 harris_env%fit_temperature
219 WRITE (unit_nr, '(T2,A,T61,ES20.8)') "Harris density fit: Entropy weight", &
220 harris_env%fit_relative_entropy_weight
221 END SELECT
222 WRITE (unit_nr, '(T2,A,T61,L20)') "Direct fitted-DM energy evaluation", &
223 harris_env%direct_density_matrix_energy
224 END SELECT
225 IF (harris_env%density_source == hden_atomic) THEN
226 WRITE (unit_nr, '(T2,A,T71,A10)') "Harris model density: Basis type", &
227 adjustr(trim(harris_env%rhoin%basis_type))
228 WRITE (unit_nr, '(T2,A,T71,I10)') "Harris model density: Number of basis functions", &
229 harris_env%rhoin%nbas
230 END IF
231 ! orbital basis
232 SELECT CASE (harris_env%orbital_basis)
233 CASE (horb_default)
234 WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model basis: ", "Atomic kind orbitals"
235 END SELECT
236
237 WRITE (unit_nr, '(T2,A)') repeat("-", 79)
238 WRITE (unit_nr, '()')
239
240 END IF ! unit_nr
241
242 CALL timestop(handle)
243
244 END SUBROUTINE harris_write_input
245
246! **************************************************************************************************
247!> \brief ...
248!> \param qs_env ...
249!> \param harris_env ...
250! **************************************************************************************************
251 SUBROUTINE harris_density_update(qs_env, harris_env)
252 TYPE(qs_environment_type), POINTER :: qs_env
253 TYPE(harris_type), POINTER :: harris_env
254
255 CHARACTER(LEN=*), PARAMETER :: routinen = 'harris_density_update'
256
257 INTEGER :: handle, i, ikind, ngto, nkind, nset, nsgf
258 INTEGER, DIMENSION(:), POINTER :: lmax, npgf
259 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: coef
260 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: density
261 REAL(kind=dp), DIMENSION(:), POINTER :: norm
262 REAL(kind=dp), DIMENSION(:, :), POINTER :: zet
263 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcc
264 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
265 TYPE(atomic_kind_type), POINTER :: atomic_kind
266 TYPE(gto_basis_set_type), POINTER :: basis_set
267 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
268 TYPE(qs_kind_type), POINTER :: qs_kind
269
270 CALL timeset(routinen, handle)
271
272 SELECT CASE (harris_env%density_source)
273 CASE (hden_atomic)
274 IF (.NOT. harris_env%rhoin%frozen) THEN
275 CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
276 nkind=nkind)
277 DO ikind = 1, nkind
278 atomic_kind => atomic_kind_set(ikind)
279 qs_kind => qs_kind_set(ikind)
280 CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, &
281 basis_type=harris_env%rhoin%basis_type)
282 CALL get_gto_basis_set(gto_basis_set=basis_set, nset=nset, lmax=lmax, nsgf=nsgf, &
283 npgf=npgf, norm_cgf=norm, zet=zet, gcc=gcc)
284 IF (nset /= 1 .OR. lmax(1) /= 0 .OR. npgf(1) /= nsgf) THEN
285 cpabort("RHOIN illegal basis type")
286 END IF
287 DO i = 1, npgf(1)
288 IF (sum(abs(gcc(1:npgf(1), i, 1))) /= maxval(abs(gcc(1:npgf(1), i, 1)))) THEN
289 cpabort("RHOIN illegal basis type")
290 END IF
291 END DO
292 !
293 ngto = npgf(1)
294 ALLOCATE (density(ngto, 2))
295 density(1:ngto, 1) = zet(1:ngto, 1)
296 density(1:ngto, 2) = 0.0_dp
297 CALL calculate_atomic_density(density, atomic_kind, qs_kind, ngto, &
298 optbasis=.false., confine=.true.)
299 ALLOCATE (coef(ngto))
300 DO i = 1, ngto
301 coef(i) = density(i, 2)/gcc(i, i, 1)/norm(i)
302 END DO
303 IF (harris_env%rhoin%nspin == 2) THEN
304 DO i = 1, SIZE(harris_env%rhoin%rhovec(ikind, 1)%rvecs, 2)
305 harris_env%rhoin%rhovec(ikind, 1)%rvecs(1:ngto, i) = coef(1:ngto)*0.5_dp
306 harris_env%rhoin%rhovec(ikind, 2)%rvecs(1:ngto, i) = coef(1:ngto)*0.5_dp
307 END DO
308 ELSE
309 DO i = 1, SIZE(harris_env%rhoin%rhovec(ikind, 1)%rvecs, 2)
310 harris_env%rhoin%rhovec(ikind, 1)%rvecs(1:ngto, i) = coef(1:ngto)
311 END DO
312 END IF
313 DEALLOCATE (density, coef)
314 END DO
315 harris_env%rhoin%frozen = .true.
316 END IF
317 CASE (hden_cube, hden_cube_fit)
318 IF (harris_env%rhoin%nspin /= 1) THEN
319 cpabort("Harris cube densities currently require a spin-restricted calculation")
320 END IF
321 IF (len_trim(harris_env%density_filename) == 0) THEN
322 cpabort("HARRIS_METHOD%FILE_DENSITY is required for cube density sources")
323 END IF
324 IF (harris_env%density_source == hden_cube_fit) THEN
325 IF (harris_env%fit_max_iter < 1) cpabort("HARRIS_METHOD%FIT_MAX_ITER has to be positive")
326 IF (harris_env%fit_eps <= 0.0_dp) cpabort("HARRIS_METHOD%FIT_EPS has to be positive")
327 IF (harris_env%fit_step_size <= 0.0_dp) THEN
328 cpabort("HARRIS_METHOD%FIT_STEP_SIZE has to be positive")
329 END IF
330 IF (harris_env%fit_max_backtrack < 0) THEN
331 cpabort("HARRIS_METHOD%FIT_MAX_BACKTRACK cannot be negative")
332 END IF
333 IF (harris_env%fit_method == hfit_relative_entropy) THEN
334 IF (harris_env%fit_temperature <= 0.0_dp) THEN
335 cpabort("HARRIS_METHOD%FIT_TEMPERATURE has to be positive for RELATIVE_ENTROPY")
336 END IF
337 IF (harris_env%fit_relative_entropy_weight < 0.0_dp) THEN
338 cpabort("HARRIS_METHOD%FIT_RELATIVE_ENTROPY_WEIGHT cannot be negative")
339 END IF
340 END IF
341 END IF
342 CASE DEFAULT
343 cpabort("Illegal value of harris_env%density_source")
344 END SELECT
345 IF (harris_env%direct_density_matrix_energy .AND. &
346 harris_env%density_source /= hden_cube_fit) THEN
347 cpabort("HARRIS_METHOD%DIRECT_DENSITY_MATRIX_ENERGY requires DENSITY_SOURCE CUBE_FIT")
348 END IF
349
350 CALL timestop(handle)
351
352 END SUBROUTINE harris_density_update
353
354! **************************************************************************************************
355!> \brief ...
356!> \param qs_env ...
357!> \param harris_env ...
358!> \param rho_struct ...
359! **************************************************************************************************
360 SUBROUTINE calculate_harris_density(qs_env, harris_env, rho_struct)
361 TYPE(qs_environment_type), POINTER :: qs_env
362 TYPE(harris_type), POINTER :: harris_env
363 TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
364
365 REAL(kind=dp), DIMENSION(:), POINTER :: tot_rho_r
366 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_gspace
367 TYPE(pw_grid_type), POINTER :: pw_grid
368 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_rspace
369
370 NULLIFY (pw_grid, rho_gspace, rho_rspace, tot_rho_r)
371
372 SELECT CASE (harris_env%density_source)
373 CASE (hden_atomic)
374 CALL calculate_harris_atomic_density(qs_env, harris_env%rhoin, rho_struct)
375 CASE (hden_cube)
376 IF (harris_env%rhoin%nspin /= 1) THEN
377 cpabort("Harris cube densities currently require a spin-restricted calculation")
378 END IF
379 IF (len_trim(harris_env%density_filename) == 0) THEN
380 cpabort("HARRIS_METHOD%FILE_DENSITY is required for DENSITY_SOURCE CUBE")
381 END IF
382 CALL read_cube_density(qs_env, rho_struct, trim(harris_env%density_filename), &
383 total_density_sign=-1, source_label="HARRIS")
384 CASE (hden_cube_fit)
385 IF (harris_env%fit_method == hfit_relative_entropy) THEN
386 IF (.NOT. harris_env%density_target_ready) THEN
387 CALL read_cube_density(qs_env, rho_struct, trim(harris_env%density_filename), &
388 total_density_sign=-1, source_label="HARRIS")
389 CALL qs_rho_get(rho_struct, rho_r=rho_rspace)
390 pw_grid => rho_rspace(1)%pw_grid
391 CALL harris_env%density_target_rspace%create(pw_grid)
392 CALL pw_copy(rho_rspace(1), harris_env%density_target_rspace)
393 harris_env%density_target_ready = .true.
394 ELSE
395 CALL qs_rho_get(rho_struct, rho_r=rho_rspace, rho_g=rho_gspace, &
396 tot_rho_r=tot_rho_r)
397 IF (harris_env%density_fit_ready) THEN
398 CALL pw_copy(harris_env%density_fit_rspace, rho_rspace(1))
399 ELSE
400 CALL pw_copy(harris_env%density_target_rspace, rho_rspace(1))
401 END IF
402 CALL pw_transfer(rho_rspace(1), rho_gspace(1))
403 tot_rho_r(1) = pw_integrate_function(rho_rspace(1), isign=-1)
404 END IF
405 ELSE IF (.NOT. harris_env%density_fit_ready) THEN
406 CALL read_cube_density(qs_env, rho_struct, trim(harris_env%density_filename), &
407 total_density_sign=-1, source_label="HARRIS")
408 IF (harris_env%direct_density_matrix_energy) THEN
409 CALL qs_rho_get(rho_struct, rho_r=rho_rspace)
410 pw_grid => rho_rspace(1)%pw_grid
411 CALL harris_env%density_target_rspace%create(pw_grid)
412 CALL pw_copy(rho_rspace(1), harris_env%density_target_rspace)
413 END IF
414 CALL fit_constrained_density(qs_env, rho_struct, harris_env%fit_max_iter, &
415 harris_env%fit_eps, harris_env%fit_step_size, &
416 harris_env%fit_max_backtrack)
417 CALL qs_rho_get(rho_struct, rho_r=rho_rspace)
418 pw_grid => rho_rspace(1)%pw_grid
419 CALL harris_env%density_fit_rspace%create(pw_grid)
420 CALL pw_copy(rho_rspace(1), harris_env%density_fit_rspace)
421 harris_env%density_fit_ready = .true.
422 ELSE
423 CALL qs_rho_get(rho_struct, rho_r=rho_rspace, rho_g=rho_gspace, &
424 tot_rho_r=tot_rho_r)
425 CALL pw_copy(harris_env%density_fit_rspace, rho_rspace(1))
426 CALL pw_transfer(rho_rspace(1), rho_gspace(1))
427 tot_rho_r(1) = pw_integrate_function(rho_rspace(1), isign=-1)
428 END IF
429 CASE DEFAULT
430 cpabort("Illegal value of harris_env%density_source")
431 END SELECT
432
433 END SUBROUTINE calculate_harris_density
434
435! **************************************************************************************************
436!> \brief Collocates an atom-centered Harris input density
437!> \param qs_env ...
438!> \param rhoin ...
439!> \param rho_struct ...
440! **************************************************************************************************
441 SUBROUTINE calculate_harris_atomic_density(qs_env, rhoin, rho_struct)
442 TYPE(qs_environment_type), POINTER :: qs_env
443 TYPE(harris_rhoin_type), INTENT(IN) :: rhoin
444 TYPE(qs_rho_type), INTENT(INOUT) :: rho_struct
445
446 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_harris_atomic_density'
447
448 INTEGER :: handle, i1, i2, iatom, ikind, ilocal, &
449 ispin, n, nkind, nlocal, nspin
450 REAL(kind=dp) :: eps_rho_rspace
451 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: vector
452 REAL(kind=dp), DIMENSION(:), POINTER :: total_rho
453 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
454 TYPE(cell_type), POINTER :: cell
455 TYPE(dft_control_type), POINTER :: dft_control
456 TYPE(distribution_1d_type), POINTER :: local_particles
457 TYPE(mp_para_env_type), POINTER :: para_env
458 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
459 TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_gspace
460 TYPE(pw_env_type), POINTER :: pw_env
461 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_rspace
462 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
463
464 CALL timeset(routinen, handle)
465
466 CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
467 eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
468 CALL get_qs_env(qs_env, &
469 atomic_kind_set=atomic_kind_set, particle_set=particle_set, &
470 local_particles=local_particles, &
471 qs_kind_set=qs_kind_set, cell=cell, pw_env=pw_env)
472
473 CALL qs_rho_get(rho_struct, rho_r=rho_rspace, rho_g=rho_gspace, &
474 tot_rho_r=total_rho)
475
476 ALLOCATE (vector(rhoin%nbas))
477
478 nkind = SIZE(rhoin%rhovec, 1)
479 nspin = SIZE(rhoin%rhovec, 2)
480
481 DO ispin = 1, nspin
482 vector = 0.0_dp
483 DO ikind = 1, nkind
484 nlocal = local_particles%n_el(ikind)
485 DO ilocal = 1, nlocal
486 iatom = local_particles%list(ikind)%array(ilocal)
487 i1 = rhoin%basptr(iatom, 1)
488 i2 = rhoin%basptr(iatom, 2)
489 n = i2 - i1 + 1
490 vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal)
491 END DO
492 END DO
493 CALL para_env%sum(vector)
494 !
495 CALL collocate_function(vector, rho_rspace(ispin), rho_gspace(ispin), &
496 atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, &
497 eps_rho_rspace, rhoin%basis_type)
498 total_rho(ispin) = pw_integrate_function(rho_rspace(ispin), isign=-1)
499 END DO
500
501 DEALLOCATE (vector)
502
503 CALL timestop(handle)
504
505 END SUBROUTINE calculate_harris_atomic_density
506
507! **************************************************************************************************
508!> \brief ...
509!> \param qs_env ...
510!> \param rhoin ...
511!> \param v_rspace ...
512!> \param calculate_forces ...
513! **************************************************************************************************
514 SUBROUTINE calculate_harris_integrals(qs_env, rhoin, v_rspace, calculate_forces)
515 TYPE(qs_environment_type), POINTER :: qs_env
516 TYPE(harris_rhoin_type), INTENT(INOUT) :: rhoin
517 TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: v_rspace
518 LOGICAL, INTENT(IN) :: calculate_forces
519
520 CHARACTER(LEN=*), PARAMETER :: routinen = 'calculate_harris_integrals'
521
522 INTEGER :: handle, i1, i2, iatom, ikind, ilocal, &
523 ispin, n, nkind, nlocal, nspin
524 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: integral, vector
525 TYPE(distribution_1d_type), POINTER :: local_particles
526 TYPE(mp_para_env_type), POINTER :: para_env
527
528 CALL timeset(routinen, handle)
529
530 CALL get_qs_env(qs_env, para_env=para_env, local_particles=local_particles)
531
532 ALLOCATE (vector(rhoin%nbas))
533 ALLOCATE (integral(rhoin%nbas))
534
535 nkind = SIZE(rhoin%rhovec, 1)
536 nspin = SIZE(rhoin%rhovec, 2)
537
538 DO ispin = 1, nspin
539 vector = 0.0_dp
540 integral = 0.0_dp
541 DO ikind = 1, nkind
542 nlocal = local_particles%n_el(ikind)
543 DO ilocal = 1, nlocal
544 iatom = local_particles%list(ikind)%array(ilocal)
545 i1 = rhoin%basptr(iatom, 1)
546 i2 = rhoin%basptr(iatom, 2)
547 n = i2 - i1 + 1
548 vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal)
549 END DO
550 END DO
551 CALL para_env%sum(vector)
552 !
553 CALL integrate_function(qs_env, v_rspace(ispin), vector, integral, &
554 calculate_forces, rhoin%basis_type)
555 DO ikind = 1, nkind
556 nlocal = local_particles%n_el(ikind)
557 DO ilocal = 1, nlocal
558 iatom = local_particles%list(ikind)%array(ilocal)
559 i1 = rhoin%basptr(iatom, 1)
560 i2 = rhoin%basptr(iatom, 2)
561 n = i2 - i1 + 1
562 rhoin%intvec(ikind, ispin)%rvecs(1:n, ilocal) = integral(i1:i2)
563 END DO
564 END DO
565 END DO
566
567 DEALLOCATE (vector, integral)
568
569 CALL timestop(handle)
570
571 END SUBROUTINE calculate_harris_integrals
572
573! **************************************************************************************************
574!> \brief ...
575!> \param harris_env ...
576!> \param vh_rspace ...
577!> \param vxc_rspace ...
578! **************************************************************************************************
579 SUBROUTINE harris_set_potentials(harris_env, vh_rspace, vxc_rspace)
580 TYPE(harris_type), POINTER :: harris_env
581 TYPE(pw_r3d_rs_type), INTENT(IN) :: vh_rspace
582 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace
583
584 INTEGER :: iab, ispin, nspins
585 TYPE(pw_grid_type), POINTER :: pw_grid
586
587 ! release possible old potentials
588 IF (ASSOCIATED(harris_env%vh_rspace%pw_grid)) THEN
589 CALL harris_env%vh_rspace%release()
590 END IF
591 IF (ASSOCIATED(harris_env%vxc_rspace)) THEN
592 DO iab = 1, SIZE(harris_env%vxc_rspace)
593 CALL harris_env%vxc_rspace(iab)%release()
594 END DO
595 DEALLOCATE (harris_env%vxc_rspace)
596 END IF
597
598 ! generate new potential data structures
599 nspins = harris_env%rhoin%nspin
600 ALLOCATE (harris_env%vxc_rspace(nspins))
601
602 pw_grid => vh_rspace%pw_grid
603 CALL harris_env%vh_rspace%create(pw_grid)
604 DO ispin = 1, nspins
605 CALL harris_env%vxc_rspace(ispin)%create(pw_grid)
606 END DO
607
608 ! copy potentials
609 CALL pw_transfer(vh_rspace, harris_env%vh_rspace)
610 IF (ASSOCIATED(vxc_rspace)) THEN
611 DO ispin = 1, nspins
612 CALL pw_transfer(vxc_rspace(ispin), harris_env%vxc_rspace(ispin))
613 CALL pw_scale(harris_env%vxc_rspace(ispin), vxc_rspace(ispin)%pw_grid%dvol)
614 END DO
615 ELSE
616 DO ispin = 1, nspins
617 CALL pw_zero(harris_env%vxc_rspace(ispin))
618 END DO
619 END IF
620
621 END SUBROUTINE harris_set_potentials
622
623END MODULE qs_harris_utils
pw_methods::pw_copy
Definition pw_methods.F:146
pw_methods::pw_integrate_function
Definition pw_methods.F:116
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
atom_kind_orbitals
calculate the orbitals for a given atomic kind type
Definition atom_kind_orbitals.F:11
atom_kind_orbitals::calculate_atomic_density
subroutine, public calculate_atomic_density(density, atomic_kind, qs_kind, ngto, iunit, optbasis, allelectron, confine)
...
Definition atom_kind_orbitals.F:571
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_unit_nr
recursive integer function, public cp_logger_get_default_unit_nr(logger, local, skip_not_ionode)
asks the default unit number of the given logger. try to use cp_logger_get_unit_nr
Definition cp_log_handling.F:576
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::hden_cube
integer, parameter, public hden_cube
Definition input_constants.F:1225
input_constants::hfit_least_squares
integer, parameter, public hfit_least_squares
Definition input_constants.F:1228
input_constants::hfit_relative_entropy
integer, parameter, public hfit_relative_entropy
Definition input_constants.F:1228
input_constants::hfun_harris
integer, parameter, public hfun_harris
Definition input_constants.F:1224
input_constants::horb_default
integer, parameter, public horb_default
Definition input_constants.F:1230
input_constants::hden_cube_fit
integer, parameter, public hden_cube_fit
Definition input_constants.F:1225
input_constants::hden_atomic
integer, parameter, public hden_atomic
Definition input_constants.F:1225
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
pw_env_types
container for various plainwaves related things
Definition pw_env_types.F:14
pw_grid_types
Definition pw_grid_types.F:13
pw_methods
Definition pw_methods.F:25
pw_types
Definition pw_types.F:24
qs_collocate_density
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
Definition qs_collocate_density.F:38
qs_collocate_density::collocate_function
subroutine, public collocate_function(vector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, particle_set, pw_env, eps_rho_rspace, basis_type)
maps a given function on the grid
Definition qs_collocate_density.F:2931
qs_density_fit
Constrained fitting of a real-space density by an AO density matrix. This is a dense reference implem...
Definition qs_density_fit.F:13
qs_density_fit::fit_constrained_density
subroutine, public fit_constrained_density(qs_env, rho_struct, max_iter, eps_rms, step_size, max_backtrack)
Fits the density currently stored in rho_struct with a spin-restricted AO density matrix....
Definition qs_density_fit.F:86
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:544
qs_external_density
Routines to handle an external density The external density can be generic and is provided by user in...
Definition qs_external_density.F:12
qs_external_density::read_cube_density
subroutine, public read_cube_density(qs_env, rho_target, filename, total_density_sign, source_label)
Read an electron density from a Gaussian cube file into a QS density grid.
Definition qs_external_density.F:148
qs_harris_types
Types needed for a for a Harris model calculation.
Definition qs_harris_types.F:14
qs_harris_utils
Harris method environment setup and handling.
Definition qs_harris_utils.F:14
qs_harris_utils::harris_write_input
subroutine, public harris_write_input(harris_env)
Print out the Harris method input section.
Definition qs_harris_utils.F:168
qs_harris_utils::harris_density_update
subroutine, public harris_density_update(qs_env, harris_env)
...
Definition qs_harris_utils.F:252
qs_harris_utils::harris_set_potentials
subroutine, public harris_set_potentials(harris_env, vh_rspace, vxc_rspace)
...
Definition qs_harris_utils.F:580
qs_harris_utils::harris_env_create
subroutine, public harris_env_create(qs_env, harris_env, harris_section)
Allocates and intitializes harris_env.
Definition qs_harris_utils.F:81
qs_harris_utils::calculate_harris_density
subroutine, public calculate_harris_density(qs_env, harris_env, rho_struct)
...
Definition qs_harris_utils.F:361
qs_integrate_potential
Integrate single or product functions over a potential on a RS grid.
Definition qs_integrate_potential.F:22
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:761
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
particle_types::particle_type
Definition particle_types.F:35
pw_env_types::pw_env_type
contained for different pw related things
Definition pw_env_types.F:70
pw_grid_types::pw_grid_type
Definition pw_grid_types.F:53
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:224
qs_harris_types::harris_rhoin_type
Definition qs_harris_types.F:38
qs_harris_types::harris_type
Contains information on the Harris method.
Definition qs_harris_types.F:66
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62