Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 2408 : 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 2408 : natom = SIZE(charges, 1)
97 2408 : nspin = SIZE(charges, 2)
98 2408 : WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
99 2408 : WRITE (UNIT=unit_nr, FMT="(T28,A)") "Hirshfeld Charges"
100 2408 : IF (nspin == 1) THEN
101 : WRITE (UNIT=unit_nr, FMT="(/,T3,A,A)") &
102 2019 : "#Atom Element Kind ", " Ref Charge Population Net charge"
103 : ELSE
104 : WRITE (UNIT=unit_nr, FMT="(/,T3,A,A)") &
105 389 : "#Atom Element Kind ", " Ref Charge Population Spin moment Net charge"
106 : END IF
107 2408 : tc1 = 0.0_dp
108 12446 : DO iatom = 1, natom
109 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
110 10038 : element_symbol=element_symbol, kind_number=ikind)
111 10038 : refc = hirshfeld_env%charges(iatom)
112 10038 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
113 10038 : 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 8571 : 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 1467 : iatom, element_symbol, ikind, refc, charges(iatom, 1), charges(iatom, 2), &
119 5868 : charges(iatom, 1) - charges(iatom, 2), zeff - SUM(charges(iatom, :))
120 : END IF
121 33989 : tc1 = tc1 + (zeff - SUM(charges(iatom, :)))
122 : END DO
123 2408 : WRITE (UNIT=unit_nr, FMT="(/,T3,A,T72,F8.3)") "Total Charge ", tc1
124 2408 : WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
125 :
126 2408 : 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 4788 : 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 4788 : NULLIFY (results)
148 4788 : CALL get_qs_env(qs_env, results=results)
149 :
150 4788 : natom = SIZE(charges, 1)
151 14364 : ALLOCATE (charges_save(natom))
152 :
153 24808 : DO iatom = 1, natom
154 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
155 20020 : kind_number=ikind)
156 20020 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
157 47762 : charges_save(iatom) = zeff - SUM(charges(iatom, :))
158 : END DO
159 :
160 : ! Store charges in results
161 4788 : description = "[HIRSHFELD-CHARGES]"
162 4788 : CALL cp_results_erase(results=results, description=description)
163 : CALL put_results(results=results, description=description, &
164 4788 : values=charges_save)
165 :
166 4788 : DEALLOCATE (charges_save)
167 :
168 4788 : 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 4970 : 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 4970 : CPASSERT(ASSOCIATED(hirshfeld_env))
196 :
197 4970 : nkind = SIZE(qs_kind_set)
198 23572 : ALLOCATE (hirshfeld_env%kind_shape_fn(nkind))
199 :
200 4970 : SELECT CASE (hirshfeld_env%shape_function_type)
201 : CASE (shape_function_gaussian)
202 13578 : DO ikind = 1, nkind
203 8626 : hirshfeld_env%kind_shape_fn(ikind)%numexp = 1
204 8626 : ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%zet(1))
205 8626 : ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%coef(1))
206 8626 : CALL get_qs_kind(qs_kind_set(ikind), element_symbol=esym)
207 8626 : rco = 2.0_dp
208 8634 : 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) &
224 0 : rco = cp_unit_to_cp2k(rco, "angstrom")
225 : END IF
226 : CASE (radius_covalent)
227 8334 : CALL get_ptable_info(symbol=esym, covalent_radius=rco, found=found)
228 8334 : IF (.NOT. found) THEN
229 0 : rco = MAX(rco, 1.0_dp)
230 : ELSE
231 8334 : IF (hirshfeld_env%use_bohr) &
232 0 : rco = cp_unit_to_cp2k(rco, "angstrom")
233 : END IF
234 : CASE (radius_single)
235 4 : CPASSERT(PRESENT(radius))
236 8630 : rco = radius
237 : END SELECT
238 8626 : al = 0.5_dp/rco**2
239 8626 : hirshfeld_env%kind_shape_fn(ikind)%zet(1) = al
240 13578 : hirshfeld_env%kind_shape_fn(ikind)%coef(1) = (al/pi)**1.5_dp
241 : END DO
242 : CASE (shape_function_density)
243 : ! calculate atomic density
244 54 : DO ikind = 1, nkind
245 36 : atomic_kind => atomic_kind_set(ikind)
246 36 : qs_kind => qs_kind_set(ikind)
247 : CALL calculate_atomic_density(ppdens(:, :), atomic_kind, qs_kind, ngto, &
248 36 : confine=.FALSE.)
249 36 : hirshfeld_env%kind_shape_fn(ikind)%numexp = ngto
250 36 : ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%zet(ngto))
251 36 : ALLOCATE (hirshfeld_env%kind_shape_fn(ikind)%coef(ngto))
252 324 : hirshfeld_env%kind_shape_fn(ikind)%zet(:) = ppdens(:, 1)
253 36 : CALL get_qs_kind(qs_kind, zeff=zeff)
254 342 : hirshfeld_env%kind_shape_fn(ikind)%coef(:) = ppdens(:, 2)/zeff
255 : END DO
256 :
257 : CASE DEFAULT
258 4970 : CPABORT("Unknown shape function")
259 : END SELECT
260 :
261 4970 : END SUBROUTINE create_shape_function
262 :
263 : ! **************************************************************************************************
264 : !> \brief ...
265 : !> \param qs_env ...
266 : !> \param hirshfeld_env ...
267 : !> \param charges ...
268 : ! **************************************************************************************************
269 4766 : SUBROUTINE comp_hirshfeld_charges(qs_env, hirshfeld_env, charges)
270 : TYPE(qs_environment_type), POINTER :: qs_env
271 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
272 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges
273 :
274 : INTEGER :: is
275 : LOGICAL :: rho_r_valid
276 : REAL(KIND=dp) :: tnfun
277 : TYPE(pw_env_type), POINTER :: pw_env
278 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
279 : TYPE(pw_r3d_rs_type) :: rhonorm
280 4766 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
281 : TYPE(qs_rho_type), POINTER :: rho
282 :
283 4766 : NULLIFY (rho_r)
284 : ! normalization function on grid
285 4766 : CALL calculate_hirshfeld_normalization(qs_env, hirshfeld_env)
286 : ! check normalization
287 4766 : tnfun = pw_integrate_function(hirshfeld_env%fnorm)
288 24720 : tnfun = ABS(tnfun - SUM(hirshfeld_env%charges))
289 : !
290 4766 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho)
291 4766 : CALL qs_rho_get(rho, rho_r=rho_r, rho_r_valid=rho_r_valid)
292 4766 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
293 4766 : CALL auxbas_pw_pool%create_pw(rhonorm)
294 : ! loop over spins
295 10298 : DO is = 1, SIZE(rho_r)
296 5532 : IF (rho_r_valid) THEN
297 : CALL hfun_scale(rhonorm%array, rho_r(is)%array, &
298 5532 : hirshfeld_env%fnorm%array)
299 : ELSE
300 0 : CPABORT("We need rho in real space")
301 : END IF
302 5532 : CALL hirshfeld_integration(qs_env, hirshfeld_env, rhonorm, charges(:, is))
303 33150 : charges(:, is) = charges(:, is)*hirshfeld_env%charges(:)
304 : END DO
305 4766 : CALL auxbas_pw_pool%give_back_pw(rhonorm)
306 :
307 4766 : END SUBROUTINE comp_hirshfeld_charges
308 : ! **************************************************************************************************
309 : !> \brief Calculate fout = fun1/fun2
310 : !> \param fout ...
311 : !> \param fun1 ...
312 : !> \param fun2 ...
313 : ! **************************************************************************************************
314 5734 : SUBROUTINE hfun_scale(fout, fun1, fun2)
315 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: fout
316 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: fun1, fun2
317 :
318 : REAL(KIND=dp), PARAMETER :: small = 1.0e-12_dp
319 :
320 : INTEGER :: i1, i2, i3, n1, n2, n3
321 :
322 5734 : n1 = SIZE(fout, 1)
323 5734 : n2 = SIZE(fout, 2)
324 5734 : n3 = SIZE(fout, 3)
325 5734 : CPASSERT(n1 == SIZE(fun1, 1))
326 5734 : CPASSERT(n2 == SIZE(fun1, 2))
327 5734 : CPASSERT(n3 == SIZE(fun1, 3))
328 5734 : CPASSERT(n1 == SIZE(fun2, 1))
329 5734 : CPASSERT(n2 == SIZE(fun2, 2))
330 5734 : CPASSERT(n3 == SIZE(fun2, 3))
331 :
332 268540 : DO i3 = 1, n3
333 12710376 : DO i2 = 1, n2
334 351242827 : DO i1 = 1, n1
335 350980021 : IF (fun2(i1, i2, i3) > small) THEN
336 117941558 : fout(i1, i2, i3) = fun1(i1, i2, i3)/fun2(i1, i2, i3)
337 : ELSE
338 220596627 : fout(i1, i2, i3) = 0.0_dp
339 : END IF
340 : END DO
341 : END DO
342 : END DO
343 :
344 5734 : END SUBROUTINE hfun_scale
345 :
346 : ! **************************************************************************************************
347 : !> \brief ...
348 : !> \param qs_env ...
349 : !> \param hirshfeld_env ...
350 : !> \param charges ...
351 : !> \param ounit ...
352 : ! **************************************************************************************************
353 22 : SUBROUTINE comp_hirshfeld_i_charges(qs_env, hirshfeld_env, charges, ounit)
354 : TYPE(qs_environment_type), POINTER :: qs_env
355 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
356 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout) :: charges
357 : INTEGER, INTENT(IN) :: ounit
358 :
359 : INTEGER, PARAMETER :: maxloop = 100
360 : REAL(KIND=dp), PARAMETER :: maxres = 1.0e-2_dp
361 :
362 : CHARACTER(len=3) :: yesno
363 : INTEGER :: iat, iloop, is, natom
364 : LOGICAL :: rho_r_valid
365 : REAL(KIND=dp) :: res, tnfun
366 : TYPE(pw_env_type), POINTER :: pw_env
367 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
368 : TYPE(pw_r3d_rs_type) :: rhonorm
369 22 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
370 : TYPE(qs_rho_type), POINTER :: rho
371 :
372 22 : NULLIFY (rho_r)
373 :
374 22 : natom = SIZE(charges, 1)
375 :
376 11 : IF (ounit > 0) WRITE (ounit, "(/,T2,A)") "Hirshfeld charge iterations: Residuals ..."
377 : !
378 22 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho=rho)
379 22 : CALL qs_rho_get(rho, rho_r=rho_r, rho_r_valid=rho_r_valid)
380 22 : CALL pw_env_get(pw_env=pw_env, auxbas_pw_pool=auxbas_pw_pool)
381 22 : CALL auxbas_pw_pool%create_pw(rhonorm)
382 : !
383 130 : DO iloop = 1, maxloop
384 :
385 : ! normalization function on grid
386 130 : CALL calculate_hirshfeld_normalization(qs_env, hirshfeld_env)
387 : ! check normalization
388 130 : tnfun = pw_integrate_function(hirshfeld_env%fnorm)
389 520 : tnfun = ABS(tnfun - SUM(hirshfeld_env%charges))
390 : ! loop over spins
391 332 : DO is = 1, SIZE(rho_r)
392 202 : IF (rho_r_valid) THEN
393 : CALL hfun_scale(rhonorm%array, rho_r(is)%array, &
394 202 : hirshfeld_env%fnorm%array)
395 : ELSE
396 0 : CPABORT("We need rho in real space")
397 : END IF
398 202 : CALL hirshfeld_integration(qs_env, hirshfeld_env, rhonorm, charges(:, is))
399 938 : charges(:, is) = charges(:, is)*hirshfeld_env%charges(:)
400 : END DO
401 : ! residual
402 130 : res = 0.0_dp
403 520 : DO iat = 1, natom
404 1126 : res = res + (SUM(charges(iat, :)) - hirshfeld_env%charges(iat))**2
405 : END DO
406 130 : res = SQRT(res/REAL(natom, KIND=dp))
407 130 : IF (ounit > 0) THEN
408 65 : yesno = "NO "
409 65 : IF (MOD(iloop, 10) == 0) yesno = "YES"
410 65 : WRITE (ounit, FMT="(F8.3)", ADVANCE=yesno) res
411 : END IF
412 : ! update
413 520 : DO iat = 1, natom
414 1126 : hirshfeld_env%charges(iat) = SUM(charges(iat, :))
415 : END DO
416 130 : IF (res < maxres) EXIT
417 :
418 : END DO
419 : !
420 22 : CALL auxbas_pw_pool%give_back_pw(rhonorm)
421 :
422 22 : END SUBROUTINE comp_hirshfeld_i_charges
423 :
424 : ! **************************************************************************************************
425 : !> \brief ...
426 : !> \param qs_env ...
427 : !> \param hirshfeld_env ...
428 : ! **************************************************************************************************
429 4896 : SUBROUTINE calculate_hirshfeld_normalization(qs_env, hirshfeld_env)
430 :
431 : TYPE(qs_environment_type), POINTER :: qs_env
432 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
433 :
434 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_hirshfeld_normalization'
435 :
436 : INTEGER :: atom_a, handle, iatom, iex, ikind, &
437 : ithread, j, natom, npme, nthread, &
438 : numexp, subpatch_pattern
439 4896 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
440 : REAL(KIND=dp) :: alpha, coef, eps_rho_rspace, radius
441 : REAL(KIND=dp), DIMENSION(3) :: ra
442 4896 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
443 4896 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
444 : TYPE(cell_type), POINTER :: cell
445 : TYPE(dft_control_type), POINTER :: dft_control
446 4896 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
447 : TYPE(pw_env_type), POINTER :: pw_env
448 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
449 : TYPE(pw_r3d_rs_type), POINTER :: fnorm
450 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
451 : TYPE(realspace_grid_type), POINTER :: rs_rho
452 :
453 4896 : CALL timeset(routineN, handle)
454 :
455 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
456 4896 : dft_control=dft_control, particle_set=particle_set, pw_env=pw_env)
457 : CALL pw_env_get(pw_env, auxbas_rs_desc=auxbas_rs_desc, auxbas_rs_grid=rs_rho, &
458 4896 : auxbas_pw_pool=auxbas_pw_pool)
459 : ! be careful in parallel nsmax is chosen with multigrid in mind!
460 4896 : CALL rs_grid_zero(rs_rho)
461 :
462 4896 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
463 4896 : ALLOCATE (pab(1, 1))
464 4896 : nthread = 1
465 4896 : ithread = 0
466 :
467 13442 : DO ikind = 1, SIZE(atomic_kind_set)
468 8546 : numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
469 8546 : IF (numexp <= 0) CYCLE
470 8546 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
471 25638 : ALLOCATE (cores(natom))
472 :
473 18268 : DO iex = 1, numexp
474 9722 : alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
475 9722 : coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
476 9722 : npme = 0
477 31830 : cores = 0
478 31830 : DO iatom = 1, natom
479 22108 : atom_a = atom_list(iatom)
480 22108 : ra(:) = pbc(particle_set(atom_a)%r, cell)
481 31830 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
482 : ! replicated realspace grid, split the atoms up between procs
483 21950 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
484 10975 : npme = npme + 1
485 10975 : cores(npme) = iatom
486 : END IF
487 : ELSE
488 158 : npme = npme + 1
489 158 : cores(npme) = iatom
490 : END IF
491 : END DO
492 29401 : DO j = 1, npme
493 11133 : iatom = cores(j)
494 11133 : atom_a = atom_list(iatom)
495 11133 : pab(1, 1) = hirshfeld_env%charges(atom_a)*coef
496 11133 : ra(:) = pbc(particle_set(atom_a)%r, cell)
497 11133 : subpatch_pattern = 0
498 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
499 : ra=ra, rb=ra, rp=ra, zetp=alpha, eps=eps_rho_rspace, &
500 : pab=pab, o1=0, o2=0, & ! without map_consistent
501 11133 : prefactor=1.0_dp, cutoff=0.0_dp)
502 :
503 : ! la_max==0 so set lmax_global to 0
504 : CALL collocate_pgf_product(0, alpha, 0, 0, 0.0_dp, 0, ra, &
505 : [0.0_dp, 0.0_dp, 0.0_dp], 1.0_dp, pab, 0, 0, rs_rho, &
506 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
507 20855 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
508 : END DO
509 : END DO
510 :
511 21988 : DEALLOCATE (cores)
512 : END DO
513 4896 : DEALLOCATE (pab)
514 :
515 4896 : NULLIFY (fnorm)
516 4896 : CALL get_hirshfeld_info(hirshfeld_env, fnorm=fnorm)
517 4896 : IF (ASSOCIATED(fnorm)) THEN
518 108 : CALL fnorm%release()
519 108 : DEALLOCATE (fnorm)
520 : END IF
521 4896 : ALLOCATE (fnorm)
522 4896 : CALL auxbas_pw_pool%create_pw(fnorm)
523 4896 : CALL set_hirshfeld_info(hirshfeld_env, fnorm=fnorm)
524 :
525 4896 : CALL transfer_rs2pw(rs_rho, fnorm)
526 :
527 4896 : CALL timestop(handle)
528 :
529 4896 : END SUBROUTINE calculate_hirshfeld_normalization
530 :
531 : ! **************************************************************************************************
532 : !> \brief ...
533 : !> \param qs_env ...
534 : !> \param hirshfeld_env ...
535 : !> \param rfun ...
536 : !> \param fval ...
537 : !> \param fderiv ...
538 : ! **************************************************************************************************
539 5734 : SUBROUTINE hirshfeld_integration(qs_env, hirshfeld_env, rfun, fval, fderiv)
540 :
541 : TYPE(qs_environment_type), POINTER :: qs_env
542 : TYPE(hirshfeld_type), POINTER :: hirshfeld_env
543 : TYPE(pw_r3d_rs_type) :: rfun
544 : REAL(KIND=dp), DIMENSION(:), INTENT(inout) :: fval
545 : REAL(KIND=dp), DIMENSION(:, :), INTENT(inout), &
546 : OPTIONAL :: fderiv
547 :
548 : CHARACTER(len=*), PARAMETER :: routineN = 'hirshfeld_integration'
549 :
550 : INTEGER :: atom_a, handle, iatom, iex, ikind, &
551 : ithread, j, natom, npme, nthread, &
552 : numexp
553 5734 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cores
554 5734 : INTEGER, DIMENSION(:), POINTER :: atom_list
555 : LOGICAL :: do_force
556 : REAL(KIND=dp) :: alpha, coef, dvol, eps_rho_rspace, radius
557 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
558 5734 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
559 5734 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
560 : TYPE(cell_type), POINTER :: cell
561 : TYPE(dft_control_type), POINTER :: dft_control
562 : TYPE(mp_para_env_type), POINTER :: para_env
563 5734 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
564 : TYPE(pw_env_type), POINTER :: pw_env
565 : TYPE(realspace_grid_desc_type), POINTER :: auxbas_rs_desc
566 : TYPE(realspace_grid_type), POINTER :: rs_v
567 :
568 5734 : CALL timeset(routineN, handle)
569 :
570 5734 : do_force = PRESENT(fderiv)
571 29192 : fval = 0.0_dp
572 5734 : dvol = rfun%pw_grid%dvol
573 :
574 5734 : NULLIFY (pw_env, auxbas_rs_desc)
575 5734 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
576 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_desc=auxbas_rs_desc, &
577 5734 : auxbas_rs_grid=rs_v)
578 5734 : CALL transfer_pw2rs(rs_v, rfun)
579 :
580 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, cell=cell, &
581 5734 : dft_control=dft_control, particle_set=particle_set)
582 5734 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
583 :
584 5734 : nthread = 1
585 5734 : ithread = 0
586 5734 : ALLOCATE (hab(1, 1), pab(1, 1))
587 :
588 15604 : DO ikind = 1, SIZE(atomic_kind_set)
589 9870 : numexp = hirshfeld_env%kind_shape_fn(ikind)%numexp
590 9870 : IF (numexp <= 0) CYCLE
591 9870 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
592 29610 : ALLOCATE (cores(natom))
593 :
594 21672 : DO iex = 1, numexp
595 11802 : alpha = hirshfeld_env%kind_shape_fn(ikind)%zet(iex)
596 11802 : coef = hirshfeld_env%kind_shape_fn(ikind)%coef(iex)
597 11802 : npme = 0
598 38158 : cores = 0
599 38158 : DO iatom = 1, natom
600 26356 : atom_a = atom_list(iatom)
601 26356 : ra(:) = pbc(particle_set(atom_a)%r, cell)
602 38158 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
603 : ! replicated realspace grid, split the atoms up between procs
604 26198 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
605 13099 : npme = npme + 1
606 13099 : cores(npme) = iatom
607 : END IF
608 : ELSE
609 158 : npme = npme + 1
610 158 : cores(npme) = iatom
611 : END IF
612 : END DO
613 :
614 34929 : DO j = 1, npme
615 13257 : iatom = cores(j)
616 13257 : atom_a = atom_list(iatom)
617 13257 : ra(:) = pbc(particle_set(atom_a)%r, cell)
618 13257 : pab(1, 1) = coef
619 13257 : hab(1, 1) = 0.0_dp
620 13257 : force_a(:) = 0.0_dp
621 13257 : force_b(:) = 0.0_dp
622 :
623 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
624 : ra=ra, rb=ra, rp=ra, &
625 : zetp=alpha, eps=eps_rho_rspace, &
626 : pab=pab, o1=0, o2=0, & ! without map_consistent
627 13257 : prefactor=1.0_dp, cutoff=1.0_dp)
628 :
629 : CALL integrate_pgf_product(0, alpha, 0, &
630 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
631 : rs_v, hab, pab=pab, o1=0, o2=0, &
632 : radius=radius, calculate_forces=do_force, &
633 : force_a=force_a, force_b=force_b, use_virial=.FALSE., &
634 13257 : use_subpatch=.TRUE., subpatch_pattern=0)
635 13257 : fval(atom_a) = fval(atom_a) + hab(1, 1)*dvol*coef
636 25059 : IF (do_force) THEN
637 0 : fderiv(:, atom_a) = fderiv(:, atom_a) + force_a(:)*dvol
638 : END IF
639 : END DO
640 :
641 : END DO
642 25474 : DEALLOCATE (cores)
643 :
644 : END DO
645 :
646 5734 : DEALLOCATE (hab, pab)
647 :
648 5734 : CALL get_qs_env(qs_env=qs_env, para_env=para_env)
649 52650 : CALL para_env%sum(fval)
650 :
651 5734 : CALL timestop(handle)
652 :
653 11468 : END SUBROUTINE hirshfeld_integration
654 :
655 : END MODULE hirshfeld_methods
|
