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 Calculate Hirshfeld charges and related functions
10 : !> \par History
11 : !> 11.2014 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE hirshfeld_methods
15 : USE ao_util, ONLY: exp_radius_very_extended
16 : USE atom_kind_orbitals, ONLY: calculate_atomic_density
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind
19 : USE cell_types, ONLY: cell_type,&
20 : pbc
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_result_methods, ONLY: cp_results_erase,&
23 : put_results
24 : USE cp_result_types, ONLY: cp_result_type
25 : USE cp_units, ONLY: cp_unit_to_cp2k
26 : USE grid_api, ONLY: GRID_FUNC_AB,&
27 : collocate_pgf_product,&
28 : integrate_pgf_product
29 : USE hirshfeld_types, ONLY: get_hirshfeld_info,&
30 : hirshfeld_type,&
31 : set_hirshfeld_info
32 : USE input_constants, ONLY: radius_covalent,&
33 : radius_default,&
34 : radius_single,&
35 : radius_user,&
36 : radius_vdw,&
37 : shape_function_density,&
38 : shape_function_gaussian
39 : USE kinds, ONLY: default_string_length,&
40 : dp
41 : USE mathconstants, ONLY: pi
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE particle_types, ONLY: particle_type
44 : USE periodic_table, ONLY: get_ptable_info
45 : USE pw_env_types, ONLY: pw_env_get,&
46 : pw_env_type
47 : USE pw_methods, ONLY: pw_integrate_function
48 : USE pw_pool_types, ONLY: pw_pool_type
49 : USE pw_types, ONLY: pw_r3d_rs_type
50 : USE qs_environment_types, ONLY: get_qs_env,&
51 : qs_environment_type
52 : USE qs_kind_types, ONLY: get_qs_kind,&
53 : qs_kind_type
54 : USE qs_rho_types, ONLY: qs_rho_get,&
55 : qs_rho_type
56 : USE realspace_grid_types, ONLY: realspace_grid_desc_type,&
57 : realspace_grid_type,&
58 : rs_grid_zero,&
59 : transfer_pw2rs,&
60 : transfer_rs2pw
61 : #include "./base/base_uses.f90"
62 :
63 : IMPLICIT NONE
64 : PRIVATE
65 :
66 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'hirshfeld_methods'
67 :
68 : PUBLIC :: create_shape_function, comp_hirshfeld_charges, &
69 : comp_hirshfeld_i_charges, write_hirshfeld_charges, &
70 : save_hirshfeld_charges
71 :
72 : ! **************************************************************************************************
73 :
74 : CONTAINS
75 :
76 : ! **************************************************************************************************
77 : !> \brief ...
78 : !> \param charges ...
79 : !> \param hirshfeld_env ...
80 : !> \param particle_set ...
81 : !> \param qs_kind_set ...
82 : !> \param unit_nr ...
83 : ! **************************************************************************************************
84 2667 : SUBROUTINE write_hirshfeld_charges(charges, hirshfeld_env, particle_set, &
85 : qs_kind_set, unit_nr)
86 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges
87 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
88 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
89 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
90 : INTEGER, INTENT(IN) :: unit_nr
91 :
92 : CHARACTER(len=2) :: element_symbol
93 : INTEGER :: iatom, ikind, natom, nspin
94 : REAL(KIND=dp) :: refc, tc1, zeff
95 :
96 2667 : natom = SIZE(charges, 1)
97 2667 : nspin = SIZE(charges, 2)
98 2667 : WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
99 2667 : WRITE (UNIT=unit_nr, FMT="(T28,A)") "Hirshfeld Charges"
100 2667 : IF (nspin == 1) THEN
101 : WRITE (UNIT=unit_nr, FMT="(/,T3,A,A)") &
102 2269 : "#Atom Element Kind ", " Ref Charge Population Net charge"
103 : ELSE
104 : WRITE (UNIT=unit_nr, FMT="(/,T3,A,A)") &
105 398 : "#Atom Element Kind ", " Ref Charge Population Spin moment Net charge"
106 : END IF
107 2667 : tc1 = 0.0_dp
108 13871 : DO iatom = 1, natom
109 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
110 11204 : element_symbol=element_symbol, kind_number=ikind)
111 11204 : refc = hirshfeld_env%charges(iatom)
112 11204 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
113 11204 : IF (nspin == 1) THEN
114 : WRITE (UNIT=unit_nr, FMT="(i7,T15,A2,T20,i3,T27,F8.3,T42,F8.3,T72,F8.3)") &
115 9712 : iatom, element_symbol, ikind, refc, charges(iatom, 1), zeff - charges(iatom, 1)
116 : ELSE
117 : WRITE (UNIT=unit_nr, FMT="(i7,T15,A2,T20,i3,T27,F8.3,T36,2F8.3,T61,F8.3,T72,F8.3)") &
118 1492 : iatom, element_symbol, ikind, refc, charges(iatom, 1), charges(iatom, 2), &
119 5968 : charges(iatom, 1) - charges(iatom, 2), zeff - SUM(charges(iatom, :))
120 : END IF
121 37771 : tc1 = tc1 + (zeff - SUM(charges(iatom, :)))
122 : END DO
123 2667 : WRITE (UNIT=unit_nr, FMT="(/,T3,A,T72,F8.3)") "Total Charge ", tc1
124 2667 : WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
125 :
126 2667 : END SUBROUTINE write_hirshfeld_charges
127 :
128 : ! **************************************************************************************************
129 : !> \brief saves the Hirshfeld charges to the results structure
130 : !> \param charges the calculated Hirshfeld charges
131 : !> \param particle_set the particle set
132 : !> \param qs_kind_set the kind set
133 : !> \param qs_env the environment
134 : ! **************************************************************************************************
135 5304 : SUBROUTINE save_hirshfeld_charges(charges, particle_set, qs_kind_set, qs_env)
136 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges
137 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
138 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
139 : TYPE(qs_environment_type), POINTER :: qs_env
140 :
141 : CHARACTER(LEN=default_string_length) :: description
142 : INTEGER :: iatom, ikind, natom
143 : REAL(KIND=dp) :: zeff
144 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges_save
145 : TYPE(cp_result_type), POINTER :: results
146 :
147 5304 : NULLIFY (results)
148 5304 : CALL get_qs_env(qs_env, results=results)
149 :
150 5304 : natom = SIZE(charges, 1)
151 15912 : ALLOCATE (charges_save(natom))
152 :
153 27652 : DO iatom = 1, natom
154 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
155 22348 : kind_number=ikind)
156 22348 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
157 52984 : charges_save(iatom) = zeff - SUM(charges(iatom, :))
158 : END DO
159 :
160 : ! Store charges in results
161 5304 : description = "[HIRSHFELD-CHARGES]"
162 5304 : CALL cp_results_erase(results=results, description=description)
163 : CALL put_results(results=results, description=description, &
164 5304 : values=charges_save)
165 :
166 5304 : DEALLOCATE (charges_save)
167 :
168 5304 : END SUBROUTINE save_hirshfeld_charges
169 :
170 : ! **************************************************************************************************
171 : !> \brief creates kind specific shape functions for Hirshfeld charges
172 : !> \param hirshfeld_env the env that holds information about Hirshfeld
173 : !> \param qs_kind_set the qs_kind_set
174 : !> \param atomic_kind_set the atomic_kind_set
175 : !> \param radius optional radius parameter to use for all atomic kinds
176 : !> \param radii_list optional list of radii to use for different atomic kinds
177 : ! **************************************************************************************************
178 5486 : SUBROUTINE create_shape_function(hirshfeld_env, qs_kind_set, atomic_kind_set, radius, radii_list)
179 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
180 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
181 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
182 : REAL(KIND=dp), OPTIONAL :: radius
183 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii_list
184 :
185 : INTEGER, PARAMETER :: ngto = 8
186 :
187 : CHARACTER(len=2) :: esym
188 : INTEGER :: ikind, nkind
189 : LOGICAL :: found
190 : REAL(KIND=dp) :: al, rco, zeff
191 : REAL(KIND=dp), DIMENSION(ngto, 2) :: ppdens
192 : TYPE(atomic_kind_type), POINTER :: atomic_kind
193 : TYPE(qs_kind_type), POINTER :: qs_kind
194 :
195 5486 : CPASSERT(ASSOCIATED(hirshfeld_env))
196 :
197 5486 : nkind = SIZE(qs_kind_set)
198 26050 : ALLOCATE (hirshfeld_env%kind_shape_fn(nkind))
199 :
200 5486 : SELECT CASE (hirshfeld_env%shape_function_type)
201 : CASE (shape_function_gaussian)
202 15024 : DO ikind = 1, nkind
203 9556 : hirshfeld_env%kind_shape_fn(ikind)%numexp = 1
204 9556 : ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%zet(1))
205 9556 : ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%coef(1))
206 9556 : CALL get_qs_kind(qs_kind_set(ikind), element_symbol=esym)
207 9556 : rco = 2.0_dp
208 9564 : SELECT CASE (hirshfeld_env%radius_type)
209 : CASE (radius_default)
210 8 : CALL get_ptable_info(symbol=esym, covalent_radius=rco, found=found)
211 8 : rco = MAX(rco, 1.0_dp)
212 : CASE (radius_user)
213 4 : CPASSERT(PRESENT(radii_list))
214 4 : CPASSERT(ASSOCIATED(radii_list))
215 4 : CPASSERT(SIZE(radii_list) == nkind)
216 : ! Note we assume that radii_list is correctly ordered
217 4 : rco = radii_list(ikind)
218 : CASE (radius_vdw)
219 276 : CALL get_ptable_info(symbol=esym, vdw_radius=rco, found=found)
220 276 : IF (.NOT. found) THEN
221 0 : rco = MAX(rco, 1.0_dp)
222 : ELSE
223 276 : IF (hirshfeld_env%use_bohr) THEN
224 0 : rco = cp_unit_to_cp2k(rco, "angstrom")
225 : END IF
226 : END IF
227 : CASE (radius_covalent)
228 9264 : CALL get_ptable_info(symbol=esym, covalent_radius=rco, found=found)
229 9264 : IF (.NOT. found) THEN
230 0 : rco = MAX(rco, 1.0_dp)
231 : ELSE
232 9264 : IF (hirshfeld_env%use_bohr) THEN
233 0 : rco = cp_unit_to_cp2k(rco, "angstrom")
234 : END IF
235 : END IF
236 : CASE (radius_single)
237 4 : CPASSERT(PRESENT(radius))
238 9560 : rco = radius
239 : END SELECT
240 9556 : al = 0.5_dp/rco**2
241 9556 : hirshfeld_env%kind_shape_fn(ikind)%zet(1) = al
242 15024 : hirshfeld_env%kind_shape_fn(ikind)%coef(1) = (al/pi)**1.5_dp
243 : END DO
244 : CASE (shape_function_density)
245 : ! calculate atomic density
246 54 : DO ikind = 1, nkind
247 36 : atomic_kind => atomic_kind_set(ikind)
248 36 : qs_kind => qs_kind_set(ikind)
249 : CALL calculate_atomic_density(ppdens(:, :), atomic_kind, qs_kind, ngto, &
250 36 : confine=.FALSE.)
251 36 : hirshfeld_env%kind_shape_fn(ikind)%numexp = ngto
252 36 : ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%zet(ngto))
253 36 : ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%coef(ngto))
254 324 : hirshfeld_env%kind_shape_fn(ikind)%zet(:) = ppdens(:, 1)
255 36 : CALL get_qs_kind(qs_kind, zeff=zeff)
256 342 : hirshfeld_env%kind_shape_fn(ikind)%coef(:) = ppdens(:, 2)/zeff
257 : END DO
258 :
259 : CASE DEFAULT
260 5486 : CPABORT("Unknown shape function")
261 : END SELECT
262 :
263 5486 : END SUBROUTINE create_shape_function
264 :
265 : ! **************************************************************************************************
266 : !> \brief ...
267 : !> \param qs_env ...
268 : !> \param hirshfeld_env ...
269 : !> \param charges ...
270 : ! **************************************************************************************************
271 5282 : SUBROUTINE comp_hirshfeld_charges(qs_env, hirshfeld_env, charges)
272 : TYPE(qs_environment_type), POINTER :: qs_env
273 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
274 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges
275 :
276 : INTEGER :: is
277 : LOGICAL :: rho_r_valid
278 : REAL(KIND=dp) :: tnfun
279 : TYPE(pw_env_type), POINTER :: pw_env
280 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
281 : TYPE(pw_r3d_rs_type) :: rhonorm
282 5282 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
283 : TYPE(qs_rho_type), POINTER :: rho
284 :
285 5282 : NULLIFY (rho_r)
286 : ! normalization function on grid
287 5282 : CALL calculate_hirshfeld_normalization(qs_env, hirshfeld_env)
288 : ! check normalization
289 5282 : tnfun = pw_integrate_function(hirshfeld_env%fnorm)
290 27564 : tnfun = ABS(tnfun - SUM(hirshfeld_env%charges))
291 : !
292 5282 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho)
293 5282 : CALL qs_rho_get(rho, rho_r=rho_r, rho_r_valid=rho_r_valid)
294 5282 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
295 5282 : CALL auxbas_pw_pool%create_pw(rhonorm)
296 : ! loop over spins
297 11348 : DO is = 1, SIZE(rho_r)
298 6066 : IF (rho_r_valid) THEN
299 : CALL hfun_scale(rhonorm%array, rho_r(is)%array, &
300 6066 : hirshfeld_env%fnorm%array)
301 : ELSE
302 0 : CPABORT("We need rho in real space")
303 : END IF
304 6066 : CALL hirshfeld_integration(qs_env, hirshfeld_env, rhonorm, charges(:, is))
305 36578 : charges(:, is) = charges(:, is)*hirshfeld_env%charges(:)
306 : END DO
307 5282 : CALL auxbas_pw_pool%give_back_pw(rhonorm)
308 :
309 5282 : END SUBROUTINE comp_hirshfeld_charges
310 : ! **************************************************************************************************
311 : !> \brief Calculate fout = fun1/fun2
312 : !> \param fout ...
313 : !> \param fun1 ...
314 : !> \param fun2 ...
315 : ! **************************************************************************************************
316 6268 : SUBROUTINE hfun_scale(fout, fun1, fun2)
317 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: fout
318 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: fun1, fun2
319 :
320 : REAL(KIND=dp), PARAMETER :: small = 1.0e-12_dp
321 :
322 : INTEGER :: i1, i2, i3, n1, n2, n3
323 :
324 6268 : n1 = SIZE(fout, 1)
325 6268 : n2 = SIZE(fout, 2)
326 6268 : n3 = SIZE(fout, 3)
327 6268 : CPASSERT(n1 == SIZE(fun1, 1))
328 6268 : CPASSERT(n2 == SIZE(fun1, 2))
329 6268 : CPASSERT(n3 == SIZE(fun1, 3))
330 6268 : CPASSERT(n1 == SIZE(fun2, 1))
331 6268 : CPASSERT(n2 == SIZE(fun2, 2))
332 6268 : CPASSERT(n3 == SIZE(fun2, 3))
333 :
334 296592 : DO i3 = 1, n3
335 14199558 : DO i2 = 1, n2
336 397297725 : DO i1 = 1, n1
337 397007401 : IF (fun2(i1, i2, i3) > small) THEN
338 139663568 : fout(i1, i2, i3) = fun1(i1, i2, i3)/fun2(i1, i2, i3)
339 : ELSE
340 243440867 : fout(i1, i2, i3) = 0.0_dp
341 : END IF
342 : END DO
343 : END DO
344 : END DO
345 :
346 6268 : END SUBROUTINE hfun_scale
347 :
348 : ! **************************************************************************************************
349 : !> \brief ...
350 : !> \param qs_env ...
351 : !> \param hirshfeld_env ...
352 : !> \param charges ...
353 : !> \param ounit ...
354 : ! **************************************************************************************************
355 22 : SUBROUTINE comp_hirshfeld_i_charges(qs_env, hirshfeld_env, charges, ounit)
356 : TYPE(qs_environment_type), POINTER :: qs_env
357 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
358 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges
359 : INTEGER, INTENT(IN) :: ounit
360 :
361 : INTEGER, PARAMETER :: maxloop = 100
362 : REAL(KIND=dp), PARAMETER :: maxres = 1.0e-2_dp
363 :
364 : CHARACTER(len=3) :: yesno
365 : INTEGER :: iat, iloop, is, natom
366 : LOGICAL :: rho_r_valid
367 : REAL(KIND=dp) :: res, tnfun
368 : TYPE(pw_env_type), POINTER :: pw_env
369 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
370 : TYPE(pw_r3d_rs_type) :: rhonorm
371 22 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
372 : TYPE(qs_rho_type), POINTER :: rho
373 :
374 22 : NULLIFY (rho_r)
375 :
376 22 : natom = SIZE(charges, 1)
377 :
378 11 : IF (ounit > 0) WRITE (ounit, "(/,T2,A)") "Hirshfeld charge iterations: Residuals ..."
379 : !
380 22 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho)
381 22 : CALL qs_rho_get(rho, rho_r=rho_r, rho_r_valid=rho_r_valid)
382 22 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
383 22 : CALL auxbas_pw_pool%create_pw(rhonorm)
384 : !
385 130 : DO iloop = 1, maxloop
386 :
387 : ! normalization function on grid
388 130 : CALL calculate_hirshfeld_normalization(qs_env, hirshfeld_env)
389 : ! check normalization
390 130 : tnfun = pw_integrate_function(hirshfeld_env%fnorm)
391 520 : tnfun = ABS(tnfun - SUM(hirshfeld_env%charges))
392 : ! loop over spins
393 332 : DO is = 1, SIZE(rho_r)
394 202 : IF (rho_r_valid) THEN
395 : CALL hfun_scale(rhonorm%array, rho_r(is)%array, &
396 202 : hirshfeld_env%fnorm%array)
397 : ELSE
398 0 : CPABORT("We need rho in real space")
399 : END IF
400 202 : CALL hirshfeld_integration(qs_env, hirshfeld_env, rhonorm, charges(:, is))
401 938 : charges(:, is) = charges(:, is)*hirshfeld_env%charges(:)
402 : END DO
403 : ! residual
404 130 : res = 0.0_dp
405 520 : DO iat = 1, natom
406 1126 : res = res + (SUM(charges(iat, :)) - hirshfeld_env%charges(iat))**2
407 : END DO
408 130 : res = SQRT(res/REAL(natom, KIND=dp))
409 130 : IF (ounit > 0) THEN
410 65 : yesno = "NO "
411 65 : IF (MOD(iloop, 10) == 0) yesno = "YES"
412 65 : WRITE (ounit, FMT="(F8.3)", ADVANCE=yesno) res
413 : END IF
414 : ! update
415 520 : DO iat = 1, natom
416 1126 : hirshfeld_env%charges(iat) = SUM(charges(iat, :))
417 : END DO
418 130 : IF (res < maxres) EXIT
419 :
420 : END DO
421 : !
422 22 : CALL auxbas_pw_pool%give_back_pw(rhonorm)
423 :
424 22 : END SUBROUTINE comp_hirshfeld_i_charges
425 :
426 : ! **************************************************************************************************
427 : !> \brief ...
428 : !> \param qs_env ...
429 : !> \param hirshfeld_env ...
430 : ! **************************************************************************************************
431 5412 : SUBROUTINE calculate_hirshfeld_normalization(qs_env, hirshfeld_env)
432 :
433 : TYPE(qs_environment_type), POINTER :: qs_env
434 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
435 :
436 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_hirshfeld_normalization'
437 :
438 : INTEGER :: atom_a, handle, iatom, iex, ikind, &
439 : ithread, j, natom, npme, nthread, &
440 : numexp, subpatch_pattern
441 5412 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
442 : REAL(KIND=dp) :: alpha, coef, eps_rho_rspace, radius
443 : REAL(KIND=dp), DIMENSION(3) :: ra
444 5412 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
445 5412 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
446 : TYPE(cell_type), POINTER :: cell
447 : TYPE(dft_control_type), POINTER :: dft_control
448 5412 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
449 : TYPE(pw_env_type), POINTER :: pw_env
450 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
451 : TYPE(pw_r3d_rs_type), POINTER :: fnorm
452 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
453 : TYPE(realspace_grid_type), POINTER :: rs_rho
454 :
455 5412 : CALL timeset(routineN, handle)
456 :
457 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
458 5412 : dft_control=dft_control, particle_set=particle_set, pw_env=pw_env)
459 : CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, auxbas_rs_grid=rs_rho, &
460 5412 : auxbas_pw_pool=auxbas_pw_pool)
461 : ! be careful in parallel nsmax is chosen with multigrid in mind!
462 5412 : CALL rs_grid_zero(rs_rho)
463 :
464 5412 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
465 5412 : ALLOCATE (pab(1, 1))
466 5412 : nthread = 1
467 5412 : ithread = 0
468 :
469 14888 : DO ikind = 1, SIZE(atomic_kind_set)
470 9476 : numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
471 9476 : IF (numexp <= 0) CYCLE
472 9476 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
473 28428 : ALLOCATE (cores(natom))
474 :
475 20128 : DO iex = 1, numexp
476 10652 : alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
477 10652 : coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
478 10652 : npme = 0
479 35088 : cores = 0
480 35088 : DO iatom = 1, natom
481 24436 : atom_a = atom_list(iatom)
482 24436 : ra(:) = pbc(particle_set(atom_a)%r, cell)
483 35088 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
484 : ! replicated realspace grid, split the atoms up between procs
485 24274 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
486 12137 : npme = npme + 1
487 12137 : cores(npme) = iatom
488 : END IF
489 : ELSE
490 162 : npme = npme + 1
491 162 : cores(npme) = iatom
492 : END IF
493 : END DO
494 32427 : DO j = 1, npme
495 12299 : iatom = cores(j)
496 12299 : atom_a = atom_list(iatom)
497 12299 : pab(1, 1) = hirshfeld_env%charges(atom_a)*coef
498 12299 : ra(:) = pbc(particle_set(atom_a)%r, cell)
499 12299 : subpatch_pattern = 0
500 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
501 : ra=ra, rb=ra, rp=ra, zetp=alpha, eps=eps_rho_rspace, &
502 : pab=pab, o1=0, o2=0, & ! without map_consistent
503 12299 : prefactor=1.0_dp, cutoff=0.0_dp)
504 :
505 : ! la_max==0 so set lmax_global to 0
506 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
507 : [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, rs_rho, &
508 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
509 22951 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
510 : END DO
511 : END DO
512 :
513 24364 : DEALLOCATE (cores)
514 : END DO
515 5412 : DEALLOCATE (pab)
516 :
517 5412 : NULLIFY (fnorm)
518 5412 : CALL get_hirshfeld_info(hirshfeld_env, fnorm=fnorm)
519 5412 : IF (ASSOCIATED(fnorm)) THEN
520 108 : CALL fnorm%release()
521 108 : DEALLOCATE (fnorm)
522 : END IF
523 5412 : ALLOCATE (fnorm)
524 5412 : CALL auxbas_pw_pool%create_pw(fnorm)
525 5412 : CALL set_hirshfeld_info(hirshfeld_env, fnorm=fnorm)
526 :
527 5412 : CALL transfer_rs2pw(rs_rho, fnorm)
528 :
529 5412 : CALL timestop(handle)
530 :
531 5412 : END SUBROUTINE calculate_hirshfeld_normalization
532 :
533 : ! **************************************************************************************************
534 : !> \brief ...
535 : !> \param qs_env ...
536 : !> \param hirshfeld_env ...
537 : !> \param rfun ...
538 : !> \param fval ...
539 : !> \param fderiv ...
540 : ! **************************************************************************************************
541 6268 : SUBROUTINE hirshfeld_integration(qs_env, hirshfeld_env, rfun, fval, fderiv)
542 :
543 : TYPE(qs_environment_type), POINTER :: qs_env
544 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
545 : TYPE(pw_r3d_rs_type) :: rfun
546 : REAL(KIND=dp), DIMENSION(:), INTENT(inout) :: fval
547 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout), &
548 : OPTIONAL :: fderiv
549 :
550 : CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_integration'
551 :
552 : INTEGER :: atom_a, handle, iatom, iex, ikind, &
553 : ithread, j, natom, npme, nthread, &
554 : numexp
555 6268 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cores
556 6268 : INTEGER, DIMENSION(:), POINTER :: atom_list
557 : LOGICAL :: do_force
558 : REAL(KIND=dp) :: alpha, coef, dvol, eps_rho_rspace, radius
559 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
560 6268 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
561 6268 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
562 : TYPE(cell_type), POINTER :: cell
563 : TYPE(dft_control_type), POINTER :: dft_control
564 : TYPE(mp_para_env_type), POINTER :: para_env
565 6268 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
566 : TYPE(pw_env_type), POINTER :: pw_env
567 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
568 : TYPE(realspace_grid_type), POINTER :: rs_v
569 :
570 6268 : CALL timeset(routineN, handle)
571 :
572 6268 : do_force = PRESENT(fderiv)
573 32104 : fval = 0.0_dp
574 6268 : dvol = rfun%pw_grid%dvol
575 :
576 6268 : NULLIFY (pw_env, auxbas_rs_desc)
577 6268 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
578 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
579 6268 : auxbas_rs_grid=rs_v)
580 6268 : CALL transfer_pw2rs(rs_v, rfun)
581 :
582 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
583 6268 : dft_control=dft_control, particle_set=particle_set)
584 6268 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
585 :
586 6268 : nthread = 1
587 6268 : ithread = 0
588 6268 : ALLOCATE (hab(1, 1), pab(1, 1))
589 :
590 17100 : DO ikind = 1, SIZE(atomic_kind_set)
591 10832 : numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
592 10832 : IF (numexp <= 0) CYCLE
593 10832 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
594 32496 : ALLOCATE (cores(natom))
595 :
596 23596 : DO iex = 1, numexp
597 12764 : alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
598 12764 : coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
599 12764 : npme = 0
600 12764 : cores = 0
601 41498 : DO iatom = 1, natom
602 28734 : atom_a = atom_list(iatom)
603 28734 : ra(:) = pbc(particle_set(atom_a)%r, cell)
604 41498 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
605 : ! replicated realspace grid, split the atoms up between procs
606 28572 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
607 14286 : npme = npme + 1
608 14286 : cores(npme) = iatom
609 : END IF
610 : ELSE
611 162 : npme = npme + 1
612 162 : cores(npme) = iatom
613 : END IF
614 : END DO
615 :
616 38044 : DO j = 1, npme
617 14448 : iatom = cores(j)
618 14448 : atom_a = atom_list(iatom)
619 14448 : ra(:) = pbc(particle_set(atom_a)%r, cell)
620 14448 : pab(1, 1) = coef
621 14448 : hab(1, 1) = 0.0_dp
622 14448 : force_a(:) = 0.0_dp
623 14448 : force_b(:) = 0.0_dp
624 :
625 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
626 : ra=ra, rb=ra, rp=ra, &
627 : zetp=alpha, eps=eps_rho_rspace, &
628 : pab=pab, o1=0, o2=0, & ! without map_consistent
629 14448 : prefactor=1.0_dp, cutoff=1.0_dp)
630 :
631 : CALL integrate_pgf_product(0, alpha, 0, &
632 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
633 : rs_v, hab, pab=pab, o1=0, o2=0, &
634 : radius=radius, calculate_forces=do_force, &
635 : force_a=force_a, force_b=force_b, use_virial=.FALSE., &
636 14448 : use_subpatch=.TRUE., subpatch_pattern=0)
637 14448 : fval(atom_a) = fval(atom_a) + hab(1, 1)*dvol*coef
638 27212 : IF (do_force) THEN
639 0 : fderiv(:, atom_a) = fderiv(:, atom_a) + force_a(:)*dvol
640 : END IF
641 : END DO
642 :
643 : END DO
644 27932 : DEALLOCATE (cores)
645 :
646 : END DO
647 :
648 6268 : DEALLOCATE (hab, pab)
649 :
650 6268 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
651 57940 : CALL para_env%sum(fval)
652 :
653 6268 : CALL timestop(handle)
654 :
655 12536 : END SUBROUTINE hirshfeld_integration
656 :
657 : END MODULE hirshfeld_methods
|
