Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Harris method 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_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 :
69 : CONTAINS
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 9032 : 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 9032 : CPASSERT(.NOT. ASSOCIATED(harris_env))
86 9032 : ALLOCATE (harris_env)
87 9032 : CALL init_harris_env(qs_env, harris_env, harris_section)
88 :
89 9032 : 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 9032 : 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 9032 : CALL timeset(routineN, handle)
111 :
112 9032 : IF (qs_env%harris_method) THEN
113 :
114 28 : CPASSERT(PRESENT(harris_section))
115 : ! get a useful output_unit
116 28 : logger => cp_get_default_logger()
117 28 : IF (logger%para_env%is_source()) THEN
118 14 : 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 28 : i_val=harris_env%energy_functional)
125 : CALL section_vals_val_get(harris_section, "DENSITY_SOURCE", &
126 28 : i_val=harris_env%density_source)
127 : CALL section_vals_val_get(harris_section, "FILE_DENSITY", &
128 28 : c_val=harris_env%density_filename)
129 : CALL section_vals_val_get(harris_section, "FIT_MAX_ITER", &
130 28 : i_val=harris_env%fit_max_iter)
131 : CALL section_vals_val_get(harris_section, "FIT_METHOD", &
132 28 : i_val=harris_env%fit_method)
133 : CALL section_vals_val_get(harris_section, "FIT_EPS", &
134 28 : r_val=harris_env%fit_eps)
135 : CALL section_vals_val_get(harris_section, "FIT_STEP_SIZE", &
136 28 : r_val=harris_env%fit_step_size)
137 : CALL section_vals_val_get(harris_section, "FIT_MAX_BACKTRACK", &
138 28 : i_val=harris_env%fit_max_backtrack)
139 : CALL section_vals_val_get(harris_section, "FIT_TEMPERATURE", &
140 28 : r_val=harris_env%fit_temperature)
141 : CALL section_vals_val_get(harris_section, "FIT_RELATIVE_ENTROPY_WEIGHT", &
142 28 : r_val=harris_env%fit_relative_entropy_weight)
143 : CALL section_vals_val_get(harris_section, "DIRECT_DENSITY_MATRIX_ENERGY", &
144 28 : l_val=harris_env%direct_density_matrix_energy)
145 : CALL section_vals_val_get(harris_section, "ORBITAL_BASIS", &
146 28 : i_val=harris_env%orbital_basis)
147 : !
148 : CALL section_vals_val_get(harris_section, "DEBUG_FORCES", &
149 28 : l_val=harris_env%debug_forces)
150 : CALL section_vals_val_get(harris_section, "DEBUG_STRESS", &
151 28 : l_val=harris_env%debug_stress)
152 :
153 : END IF
154 :
155 9032 : CALL timestop(handle)
156 :
157 9032 : 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 28 : 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 28 : CALL timeset(routineN, handle)
176 :
177 28 : logger => cp_get_default_logger()
178 28 : IF (logger%para_env%is_source()) THEN
179 14 : unit_nr = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
180 : ELSE
181 : unit_nr = -1
182 : END IF
183 :
184 14 : IF (unit_nr > 0) THEN
185 :
186 : WRITE (unit_nr, '(/,T2,A)') &
187 14 : "!"//REPEAT("-", 29)//" Harris Model "//REPEAT("-", 29)//"!"
188 :
189 : ! Type of energy functional
190 28 : SELECT CASE (harris_env%energy_functional)
191 : CASE (hfun_harris)
192 14 : WRITE (unit_nr, '(T2,A,T61,A20)') "Energy Functional: ", "Harris"
193 : END SELECT
194 : ! density source
195 18 : SELECT CASE (harris_env%density_source)
196 : CASE (hden_atomic)
197 4 : WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model density: Type", " Atomic kind density"
198 : CASE (hden_cube)
199 2 : WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model density: Type", "Cube file"
200 2 : WRITE (unit_nr, '(T2,A,T31,A)') "Harris model density: File", &
201 4 : TRIM(harris_env%density_filename)
202 : CASE (hden_cube_fit)
203 8 : WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model density: Type", "Constrained cube fit"
204 8 : WRITE (unit_nr, '(T2,A,T31,A)') "Harris model density: File", &
205 16 : TRIM(harris_env%density_filename)
206 8 : WRITE (unit_nr, '(T2,A,T61,I20)') "Harris density fit: Maximum iterations", &
207 16 : harris_env%fit_max_iter
208 8 : WRITE (unit_nr, '(T2,A,T61,ES20.8)') "Harris density fit: RMS target", &
209 16 : harris_env%fit_eps
210 8 : WRITE (unit_nr, '(T2,A,T61,ES20.8)') "Harris density fit: Initial step size", &
211 16 : harris_env%fit_step_size
212 12 : SELECT CASE (harris_env%fit_method)
213 : CASE (hfit_least_squares)
214 4 : WRITE (unit_nr, '(T2,A,T61,A20)') "Harris density fit: Objective", "Least squares"
215 : CASE (hfit_relative_entropy)
216 4 : WRITE (unit_nr, '(T2,A,T61,A20)') "Harris density fit: Objective", "Relative entropy"
217 4 : WRITE (unit_nr, '(T2,A,T61,ES20.8)') "Harris density fit: Temperature", &
218 8 : harris_env%fit_temperature
219 4 : WRITE (unit_nr, '(T2,A,T61,ES20.8)') "Harris density fit: Entropy weight", &
220 16 : harris_env%fit_relative_entropy_weight
221 : END SELECT
222 8 : WRITE (unit_nr, '(T2,A,T61,L20)') "Direct fitted-DM energy evaluation", &
223 30 : harris_env%direct_density_matrix_energy
224 : END SELECT
225 14 : IF (harris_env%density_source == hden_atomic) THEN
226 4 : WRITE (unit_nr, '(T2,A,T71,A10)') "Harris model density: Basis type", &
227 8 : ADJUSTR(TRIM(harris_env%rhoin%basis_type))
228 4 : WRITE (unit_nr, '(T2,A,T71,I10)') "Harris model density: Number of basis functions", &
229 8 : harris_env%rhoin%nbas
230 : END IF
231 : ! orbital basis
232 28 : SELECT CASE (harris_env%orbital_basis)
233 : CASE (horb_default)
234 14 : WRITE (unit_nr, '(T2,A,T61,A20)') "Harris model basis: ", "Atomic kind orbitals"
235 : END SELECT
236 :
237 14 : WRITE (unit_nr, '(T2,A)') REPEAT("-", 79)
238 14 : WRITE (unit_nr, '()')
239 :
240 : END IF ! unit_nr
241 :
242 28 : CALL timestop(handle)
243 :
244 28 : END SUBROUTINE harris_write_input
245 :
246 : ! **************************************************************************************************
247 : !> \brief ...
248 : !> \param qs_env ...
249 : !> \param harris_env ...
250 : ! **************************************************************************************************
251 72 : 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 72 : INTEGER, DIMENSION(:), POINTER :: lmax, npgf
259 72 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: coef
260 72 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: density
261 72 : REAL(KIND=dp), DIMENSION(:), POINTER :: norm
262 72 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zet
263 72 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcc
264 72 : 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 72 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
268 : TYPE(qs_kind_type), POINTER :: qs_kind
269 :
270 72 : CALL timeset(routineN, handle)
271 :
272 116 : SELECT CASE (harris_env%density_source)
273 : CASE (hden_atomic)
274 44 : 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 8 : nkind=nkind)
277 30 : DO ikind = 1, nkind
278 22 : atomic_kind => atomic_kind_set(ikind)
279 22 : qs_kind => qs_kind_set(ikind)
280 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_set, &
281 22 : 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 22 : npgf=npgf, norm_cgf=norm, zet=zet, gcc=gcc)
284 22 : IF (nset /= 1 .OR. lmax(1) /= 0 .OR. npgf(1) /= nsgf) THEN
285 0 : CPABORT("RHOIN illegal basis type")
286 : END IF
287 168 : DO i = 1, npgf(1)
288 2116 : IF (SUM(ABS(gcc(1:npgf(1), i, 1))) /= MAXVAL(ABS(gcc(1:npgf(1), i, 1)))) THEN
289 0 : CPABORT("RHOIN illegal basis type")
290 : END IF
291 : END DO
292 : !
293 22 : ngto = npgf(1)
294 66 : ALLOCATE (density(ngto, 2))
295 168 : density(1:ngto, 1) = zet(1:ngto, 1)
296 168 : density(1:ngto, 2) = 0.0_dp
297 : CALL calculate_atomic_density(density, atomic_kind, qs_kind, ngto, &
298 22 : optbasis=.FALSE., confine=.TRUE.)
299 66 : ALLOCATE (coef(ngto))
300 168 : DO i = 1, ngto
301 168 : coef(i) = density(i, 2)/gcc(i, i, 1)/norm(i)
302 : END DO
303 22 : IF (harris_env%rhoin%nspin == 2) THEN
304 10 : DO i = 1, SIZE(harris_env%rhoin%rhovec(ikind, 1)%rvecs, 2)
305 30 : harris_env%rhoin%rhovec(ikind, 1)%rvecs(1:ngto, i) = coef(1:ngto)*0.5_dp
306 36 : harris_env%rhoin%rhovec(ikind, 2)%rvecs(1:ngto, i) = coef(1:ngto)*0.5_dp
307 : END DO
308 : ELSE
309 30 : DO i = 1, SIZE(harris_env%rhoin%rhovec(ikind, 1)%rvecs, 2)
310 120 : harris_env%rhoin%rhovec(ikind, 1)%rvecs(1:ngto, i) = coef(1:ngto)
311 : END DO
312 : END IF
313 52 : DEALLOCATE (density, coef)
314 : END DO
315 8 : harris_env%rhoin%frozen = .TRUE.
316 : END IF
317 : CASE (hden_cube, hden_cube_fit)
318 28 : IF (harris_env%rhoin%nspin /= 1) THEN
319 0 : CPABORT("Harris cube densities currently require a spin-restricted calculation")
320 : END IF
321 28 : IF (LEN_TRIM(harris_env%density_filename) == 0) THEN
322 0 : CPABORT("HARRIS_METHOD%FILE_DENSITY is required for cube density sources")
323 : END IF
324 28 : IF (harris_env%density_source == hden_cube_fit) THEN
325 24 : IF (harris_env%fit_max_iter < 1) CPABORT("HARRIS_METHOD%FIT_MAX_ITER has to be positive")
326 24 : IF (harris_env%fit_eps <= 0.0_dp) CPABORT("HARRIS_METHOD%FIT_EPS has to be positive")
327 24 : IF (harris_env%fit_step_size <= 0.0_dp) THEN
328 0 : CPABORT("HARRIS_METHOD%FIT_STEP_SIZE has to be positive")
329 : END IF
330 24 : IF (harris_env%fit_max_backtrack < 0) THEN
331 0 : CPABORT("HARRIS_METHOD%FIT_MAX_BACKTRACK cannot be negative")
332 : END IF
333 24 : IF (harris_env%fit_method == hfit_relative_entropy) THEN
334 16 : IF (harris_env%fit_temperature <= 0.0_dp) THEN
335 0 : CPABORT("HARRIS_METHOD%FIT_TEMPERATURE has to be positive for RELATIVE_ENTROPY")
336 : END IF
337 16 : IF (harris_env%fit_relative_entropy_weight < 0.0_dp) THEN
338 0 : CPABORT("HARRIS_METHOD%FIT_RELATIVE_ENTROPY_WEIGHT cannot be negative")
339 : END IF
340 : END IF
341 : END IF
342 : CASE DEFAULT
343 72 : CPABORT("Illegal value of harris_env%density_source")
344 : END SELECT
345 72 : IF (harris_env%direct_density_matrix_energy .AND. &
346 : harris_env%density_source /= hden_cube_fit) THEN
347 0 : CPABORT("HARRIS_METHOD%DIRECT_DENSITY_MATRIX_ENERGY requires DENSITY_SOURCE CUBE_FIT")
348 : END IF
349 :
350 72 : CALL timestop(handle)
351 :
352 144 : END SUBROUTINE harris_density_update
353 :
354 : ! **************************************************************************************************
355 : !> \brief ...
356 : !> \param qs_env ...
357 : !> \param harris_env ...
358 : !> \param rho_struct ...
359 : ! **************************************************************************************************
360 96 : 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 96 : REAL(KIND=dp), DIMENSION(:), POINTER :: tot_rho_r
366 96 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_gspace
367 : TYPE(pw_grid_type), POINTER :: pw_grid
368 96 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_rspace
369 :
370 96 : NULLIFY (pw_grid, rho_gspace, rho_rspace, tot_rho_r)
371 :
372 164 : SELECT CASE (harris_env%density_source)
373 : CASE (hden_atomic)
374 68 : CALL calculate_harris_atomic_density(qs_env, harris_env%rhoin, rho_struct)
375 : CASE (hden_cube)
376 4 : IF (harris_env%rhoin%nspin /= 1) THEN
377 0 : CPABORT("Harris cube densities currently require a spin-restricted calculation")
378 : END IF
379 4 : IF (LEN_TRIM(harris_env%density_filename) == 0) THEN
380 0 : 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 4 : total_density_sign=-1, source_label="HARRIS")
384 : CASE (hden_cube_fit)
385 24 : IF (harris_env%fit_method == hfit_relative_entropy) THEN
386 16 : IF (.NOT. harris_env%density_target_ready) THEN
387 : CALL read_cube_density(qs_env, rho_struct, TRIM(harris_env%density_filename), &
388 8 : total_density_sign=-1, source_label="HARRIS")
389 8 : CALL qs_rho_get(rho_struct, rho_r=rho_rspace)
390 8 : pw_grid => rho_rspace(1)%pw_grid
391 8 : CALL harris_env%density_target_rspace%create(pw_grid)
392 8 : CALL pw_copy(rho_rspace(1), harris_env%density_target_rspace)
393 8 : harris_env%density_target_ready = .TRUE.
394 : ELSE
395 : CALL qs_rho_get(rho_struct, rho_r=rho_rspace, rho_g=rho_gspace, &
396 8 : tot_rho_r=tot_rho_r)
397 8 : IF (harris_env%density_fit_ready) THEN
398 8 : CALL pw_copy(harris_env%density_fit_rspace, rho_rspace(1))
399 : ELSE
400 0 : CALL pw_copy(harris_env%density_target_rspace, rho_rspace(1))
401 : END IF
402 8 : CALL pw_transfer(rho_rspace(1), rho_gspace(1))
403 8 : tot_rho_r(1) = pw_integrate_function(rho_rspace(1), isign=-1)
404 : END IF
405 8 : ELSE IF (.NOT. harris_env%density_fit_ready) THEN
406 : CALL read_cube_density(qs_env, rho_struct, TRIM(harris_env%density_filename), &
407 8 : total_density_sign=-1, source_label="HARRIS")
408 8 : IF (harris_env%direct_density_matrix_energy) THEN
409 4 : CALL qs_rho_get(rho_struct, rho_r=rho_rspace)
410 4 : pw_grid => rho_rspace(1)%pw_grid
411 4 : CALL harris_env%density_target_rspace%create(pw_grid)
412 4 : 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 8 : harris_env%fit_max_backtrack)
417 8 : CALL qs_rho_get(rho_struct, rho_r=rho_rspace)
418 8 : pw_grid => rho_rspace(1)%pw_grid
419 8 : CALL harris_env%density_fit_rspace%create(pw_grid)
420 8 : CALL pw_copy(rho_rspace(1), harris_env%density_fit_rspace)
421 8 : harris_env%density_fit_ready = .TRUE.
422 : ELSE
423 : CALL qs_rho_get(rho_struct, rho_r=rho_rspace, rho_g=rho_gspace, &
424 0 : tot_rho_r=tot_rho_r)
425 0 : CALL pw_copy(harris_env%density_fit_rspace, rho_rspace(1))
426 0 : CALL pw_transfer(rho_rspace(1), rho_gspace(1))
427 0 : tot_rho_r(1) = pw_integrate_function(rho_rspace(1), isign=-1)
428 : END IF
429 : CASE DEFAULT
430 0 : CPABORT("Illegal value of harris_env%density_source")
431 : END SELECT
432 :
433 96 : 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 68 : 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 68 : REAL(KIND=dp), DIMENSION(:), POINTER :: total_rho
453 68 : 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 68 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
459 68 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_gspace
460 : TYPE(pw_env_type), POINTER :: pw_env
461 68 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_rspace
462 68 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
463 :
464 68 : CALL timeset(routineN, handle)
465 :
466 68 : CALL get_qs_env(qs_env, dft_control=dft_control, para_env=para_env)
467 68 : 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 68 : 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 68 : tot_rho_r=total_rho)
475 :
476 204 : ALLOCATE (vector(rhoin%nbas))
477 :
478 68 : nkind = SIZE(rhoin%rhovec, 1)
479 68 : nspin = SIZE(rhoin%rhovec, 2)
480 :
481 162 : DO ispin = 1, nspin
482 94 : vector = 0.0_dp
483 374 : DO ikind = 1, nkind
484 280 : nlocal = local_particles%n_el(ikind)
485 564 : DO ilocal = 1, nlocal
486 190 : iatom = local_particles%list(ikind)%array(ilocal)
487 190 : i1 = rhoin%basptr(iatom, 1)
488 190 : i2 = rhoin%basptr(iatom, 2)
489 190 : n = i2 - i1 + 1
490 1704 : vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal)
491 : END DO
492 : END DO
493 94 : 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 94 : eps_rho_rspace, rhoin%basis_type)
498 162 : total_rho(ispin) = pw_integrate_function(rho_rspace(ispin), isign=-1)
499 : END DO
500 :
501 68 : DEALLOCATE (vector)
502 :
503 68 : CALL timestop(handle)
504 :
505 68 : 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 0 : 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 0 : CALL timeset(routineN, handle)
529 :
530 0 : CALL get_qs_env(qs_env, para_env=para_env, local_particles=local_particles)
531 :
532 0 : ALLOCATE (vector(rhoin%nbas))
533 0 : ALLOCATE (integral(rhoin%nbas))
534 :
535 0 : nkind = SIZE(rhoin%rhovec, 1)
536 0 : nspin = SIZE(rhoin%rhovec, 2)
537 :
538 0 : DO ispin = 1, nspin
539 0 : vector = 0.0_dp
540 0 : integral = 0.0_dp
541 0 : DO ikind = 1, nkind
542 0 : nlocal = local_particles%n_el(ikind)
543 0 : DO ilocal = 1, nlocal
544 0 : iatom = local_particles%list(ikind)%array(ilocal)
545 0 : i1 = rhoin%basptr(iatom, 1)
546 0 : i2 = rhoin%basptr(iatom, 2)
547 0 : n = i2 - i1 + 1
548 0 : vector(i1:i2) = rhoin%rhovec(ikind, ispin)%rvecs(1:n, ilocal)
549 : END DO
550 : END DO
551 0 : CALL para_env%sum(vector)
552 : !
553 : CALL integrate_function(qs_env, v_rspace(ispin), vector, integral, &
554 0 : calculate_forces, rhoin%basis_type)
555 0 : DO ikind = 1, nkind
556 0 : nlocal = local_particles%n_el(ikind)
557 0 : DO ilocal = 1, nlocal
558 0 : iatom = local_particles%list(ikind)%array(ilocal)
559 0 : i1 = rhoin%basptr(iatom, 1)
560 0 : i2 = rhoin%basptr(iatom, 2)
561 0 : n = i2 - i1 + 1
562 0 : rhoin%intvec(ikind, ispin)%rvecs(1:n, ilocal) = integral(i1:i2)
563 : END DO
564 : END DO
565 : END DO
566 :
567 0 : DEALLOCATE (vector, integral)
568 :
569 0 : CALL timestop(handle)
570 :
571 0 : END SUBROUTINE calculate_harris_integrals
572 :
573 : ! **************************************************************************************************
574 : !> \brief ...
575 : !> \param harris_env ...
576 : !> \param vh_rspace ...
577 : !> \param vxc_rspace ...
578 : ! **************************************************************************************************
579 116 : 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 116 : IF (ASSOCIATED(harris_env%vh_rspace%pw_grid)) THEN
589 88 : CALL harris_env%vh_rspace%release()
590 : END IF
591 116 : IF (ASSOCIATED(harris_env%vxc_rspace)) THEN
592 200 : DO iab = 1, SIZE(harris_env%vxc_rspace)
593 200 : CALL harris_env%vxc_rspace(iab)%release()
594 : END DO
595 88 : DEALLOCATE (harris_env%vxc_rspace)
596 : END IF
597 :
598 : ! generate new potential data structures
599 116 : nspins = harris_env%rhoin%nspin
600 490 : ALLOCATE (harris_env%vxc_rspace(nspins))
601 :
602 116 : pw_grid => vh_rspace%pw_grid
603 116 : CALL harris_env%vh_rspace%create(pw_grid)
604 258 : DO ispin = 1, nspins
605 258 : CALL harris_env%vxc_rspace(ispin)%create(pw_grid)
606 : END DO
607 :
608 : ! copy potentials
609 116 : CALL pw_transfer(vh_rspace, harris_env%vh_rspace)
610 116 : IF (ASSOCIATED(vxc_rspace)) THEN
611 210 : DO ispin = 1, nspins
612 118 : CALL pw_transfer(vxc_rspace(ispin), harris_env%vxc_rspace(ispin))
613 210 : CALL pw_scale(harris_env%vxc_rspace(ispin), vxc_rspace(ispin)%pw_grid%dvol)
614 : END DO
615 : ELSE
616 48 : DO ispin = 1, nspins
617 48 : CALL pw_zero(harris_env%vxc_rspace(ispin))
618 : END DO
619 : END IF
620 :
621 116 : END SUBROUTINE harris_set_potentials
622 :
623 : END MODULE qs_harris_utils
|
