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
10 : !> \author JGH (01.2026)
11 : ! **************************************************************************************************
12 : MODULE accint_weights_forces
13 : USE ao_util, ONLY: exp_radius_very_extended
14 : USE atomic_kind_types, ONLY: atomic_kind_type,&
15 : get_atomic_kind
16 : USE cell_types, ONLY: cell_type,&
17 : pbc
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
20 : USE grid_api, ONLY: GRID_FUNC_AB,&
21 : collocate_pgf_product,&
22 : integrate_pgf_product
23 : USE input_constants, ONLY: sic_none,&
24 : xc_none
25 : USE input_section_types, ONLY: section_vals_type,&
26 : section_vals_val_get
27 : USE kinds, ONLY: dp
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE orbital_pointers, ONLY: coset,&
30 : ncoset
31 : USE particle_types, ONLY: particle_type
32 : USE pw_env_types, ONLY: pw_env_get,&
33 : pw_env_type
34 : USE pw_grids, ONLY: pw_grid_compare
35 : USE pw_methods, ONLY: pw_axpy,&
36 : pw_multiply_with,&
37 : pw_scale,&
38 : pw_transfer,&
39 : pw_zero
40 : USE pw_pool_types, ONLY: pw_pool_p_type,&
41 : pw_pool_type
42 : USE pw_types, ONLY: pw_c1d_gs_type,&
43 : pw_r3d_rs_type
44 : USE qs_environment_types, ONLY: get_qs_env,&
45 : qs_environment_type
46 : USE qs_force_types, ONLY: qs_force_type
47 : USE qs_fxc, ONLY: qs_fxc_create
48 : USE qs_kind_types, ONLY: qs_kind_type
49 : USE qs_ks_types, ONLY: qs_ks_env_type
50 : USE qs_rho_atom_types, ONLY: rho_atom_type
51 : USE qs_rho_types, ONLY: qs_rho_get,&
52 : qs_rho_type
53 : USE realspace_grid_types, ONLY: realspace_grid_type,&
54 : rs_grid_zero,&
55 : transfer_pw2rs,&
56 : transfer_rs2pw
57 : USE skala_gpw_functional, ONLY: native_skala_gapw_atom_composite_requested,&
58 : native_skala_gapw_composite_reference,&
59 : native_skala_uses_atom_composite_grid,&
60 : skala_gapw_representation,&
61 : skala_gpw_weight_derivative,&
62 : xc_section_uses_native_skala_grid
63 : USE virial_types, ONLY: virial_type
64 : USE xc, ONLY: xc_exc_pw_create,&
65 : xc_vxc_pw_create
66 : USE xc_gauxc_functional, ONLY: gauxc_gapw_has_paw_pseudopotentials
67 : USE xc_input_constants, ONLY: skala_gapw_paw_one_center
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 :
74 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
75 :
76 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'accint_weights_forces'
77 :
78 : PUBLIC :: accint_weight_force, weight_function_pwgrid
79 :
80 : CONTAINS
81 :
82 : ! **************************************************************************************************
83 : !> \brief ...
84 : !> \param qs_env ...
85 : !> \param rho ...
86 : !> \param rho1 ...
87 : !> \param order ...
88 : !> \param xc_section ...
89 : !> \param triplet ...
90 : !> \param force_scale ...
91 : ! **************************************************************************************************
92 1548 : SUBROUTINE accint_weight_force(qs_env, rho, rho1, order, xc_section, triplet, force_scale)
93 : TYPE(qs_environment_type), POINTER :: qs_env
94 : TYPE(qs_rho_type), POINTER :: rho, rho1
95 : INTEGER, INTENT(IN) :: order
96 : TYPE(section_vals_type), POINTER :: xc_section
97 : LOGICAL, INTENT(IN), OPTIONAL :: triplet
98 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: force_scale
99 :
100 : CHARACTER(len=*), PARAMETER :: routineN = 'accint_weight_force'
101 :
102 : INTEGER :: atom_a, handle, i, iatom, ikind, natom, &
103 : natom_of_kind, nkind, ounit, oweight
104 1548 : INTEGER, DIMENSION(:), POINTER :: atom_list
105 : LOGICAL :: composite_reference, lr_triplet, &
106 : native_grid_diagnostics, &
107 : native_skala_grid, uf_grid, use_virial
108 : REAL(KIND=dp) :: my_force_scale
109 1548 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: calpha, cvalue
110 1548 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aforce
111 : REAL(KIND=dp), DIMENSION(3) :: tforce
112 : REAL(KIND=dp), DIMENSION(3, 3) :: avirial
113 1548 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
114 : TYPE(dft_control_type), POINTER :: dft_control
115 : TYPE(mp_para_env_type), POINTER :: para_env
116 : TYPE(pw_env_type), POINTER :: pw_env
117 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
118 : TYPE(pw_r3d_rs_type) :: e_force_rspace, e_rspace
119 1548 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
120 1548 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
121 : TYPE(qs_ks_env_type), POINTER :: ks_env
122 : TYPE(virial_type), POINTER :: virial
123 :
124 1548 : CALL timeset(routineN, handle)
125 :
126 1548 : CALL get_qs_env(qs_env, dft_control=dft_control, qs_kind_set=qs_kind_set)
127 :
128 : ! Composite references replace the PW XC quadrature and its Gaussian weight derivative.
129 : ! The public selector applies only when a pseudopotential kind actually uses PAW one-center
130 : ! data; it must not change all-electron GAPW integration.
131 : composite_reference = native_skala_gapw_composite_reference(xc_section) .OR. &
132 1548 : native_skala_gapw_atom_composite_requested(xc_section)
133 : IF (.NOT. composite_reference) THEN
134 : composite_reference = skala_gapw_representation(xc_section) == &
135 : skala_gapw_paw_one_center .AND. &
136 1548 : gauxc_gapw_has_paw_pseudopotentials(qs_kind_set)
137 : END IF
138 1520 : IF (.NOT. composite_reference .AND. xc_section_uses_native_skala_grid(xc_section)) THEN
139 24 : composite_reference = native_skala_uses_atom_composite_grid(xc_section)
140 : END IF
141 1520 : IF (composite_reference) THEN
142 52 : CALL timestop(handle)
143 52 : RETURN
144 : END IF
145 :
146 1496 : IF (dft_control%qs_control%gapw_control%accurate_xcint) THEN
147 :
148 406 : CALL get_qs_env(qs_env=qs_env, force=force, virial=virial)
149 406 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
150 :
151 406 : ounit = cp_logger_get_default_io_unit()
152 :
153 406 : CALL get_qs_env(qs_env, natom=natom, nkind=nkind)
154 1218 : ALLOCATE (aforce(3, natom))
155 1624 : ALLOCATE (calpha(nkind), cvalue(nkind))
156 1236 : cvalue = 1.0_dp
157 1236 : calpha(1:nkind) = dft_control%qs_control%gapw_control%aw(1:nkind)
158 406 : oweight = dft_control%qs_control%gapw_control%oweights
159 :
160 406 : CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env)
161 406 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, xc_pw_pool=xc_pw_pool)
162 406 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
163 406 : IF (uf_grid) THEN
164 44 : CALL xc_pw_pool%create_pw(e_rspace)
165 : ELSE
166 362 : CALL auxbas_pw_pool%create_pw(e_rspace)
167 : END IF
168 :
169 406 : lr_triplet = .FALSE.
170 406 : IF (PRESENT(triplet)) lr_triplet = triplet
171 406 : my_force_scale = 1.0_dp
172 406 : IF (PRESENT(force_scale)) my_force_scale = force_scale
173 :
174 406 : CALL xc_density(qs_env, rho, rho1, order, xc_section, lr_triplet, e_rspace)
175 :
176 406 : IF (uf_grid) THEN
177 44 : CALL auxbas_pw_pool%create_pw(e_force_rspace)
178 : BLOCK
179 : TYPE(pw_c1d_gs_type) :: e_g_aux, e_g_xc
180 44 : CALL xc_pw_pool%create_pw(e_g_xc)
181 44 : CALL auxbas_pw_pool%create_pw(e_g_aux)
182 44 : CALL pw_transfer(e_rspace, e_g_xc)
183 44 : CALL pw_transfer(e_g_xc, e_g_aux)
184 44 : CALL pw_transfer(e_g_aux, e_force_rspace)
185 44 : CALL auxbas_pw_pool%give_back_pw(e_g_aux)
186 88 : CALL xc_pw_pool%give_back_pw(e_g_xc)
187 : END BLOCK
188 44 : CALL pw_scale(e_force_rspace, e_force_rspace%pw_grid%dvol)
189 : ! integrate energy field to get force
190 44 : CALL gauss_grid_force(e_force_rspace, qs_env, oweight, calpha, cvalue, aforce, avirial)
191 : !
192 44 : CALL auxbas_pw_pool%give_back_pw(e_force_rspace)
193 : ELSE
194 362 : CALL pw_scale(e_rspace, e_rspace%pw_grid%dvol)
195 362 : CALL gauss_grid_force(e_rspace, qs_env, oweight, calpha, cvalue, aforce, avirial)
196 : END IF
197 :
198 406 : IF (uf_grid) THEN
199 44 : CALL xc_pw_pool%give_back_pw(e_rspace)
200 : ELSE
201 362 : CALL auxbas_pw_pool%give_back_pw(e_rspace)
202 : END IF
203 :
204 406 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
205 406 : native_skala_grid = xc_section_uses_native_skala_grid(xc_section)
206 406 : native_grid_diagnostics = .FALSE.
207 406 : IF (native_skala_grid) THEN
208 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%GAUXC%NATIVE_GRID_DIAGNOSTICS", &
209 0 : l_val=native_grid_diagnostics)
210 : END IF
211 1236 : DO ikind = 1, nkind
212 830 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
213 2482 : DO iatom = 1, natom_of_kind
214 1246 : atom_a = atom_list(iatom)
215 1246 : IF (native_grid_diagnostics) THEN
216 0 : IF (ounit > 0) THEN
217 : WRITE (UNIT=ounit, FMT="(T2,A,1X,I0,3(1X,ES20.12))") &
218 0 : "SKALA_GPW| Accurate-XCINT atom force", atom_a, my_force_scale*aforce(:, atom_a)
219 : END IF
220 : END IF
221 : force(ikind)%rho_elec(1:3, iatom) = &
222 5814 : force(ikind)%rho_elec(1:3, iatom) + my_force_scale*aforce(1:3, atom_a)
223 : END DO
224 : END DO
225 406 : IF (use_virial) THEN
226 754 : virial%pv_exc = virial%pv_exc + my_force_scale*avirial
227 754 : virial%pv_virial = virial%pv_virial + my_force_scale*avirial
228 : END IF
229 :
230 : IF (debug_this_module) THEN
231 : CALL get_qs_env(qs_env, para_env=para_env)
232 : CALL para_env%sum(aforce)
233 : IF (ounit > 0) THEN
234 : WRITE (UNIT=ounit, FMT="(/,T2,A)") "ACCINT| Debug Accurate-XCINT atom force [a.u.]"
235 : DO i = 1, 3
236 : tforce(i) = my_force_scale*SUM(aforce(i, :))
237 : END DO
238 : DO iatom = 1, natom
239 : WRITE (UNIT=ounit, FMT="(T2,A,1X,I0,T24,3(1X,F18.12))") &
240 : "ACCINT| ", iatom, my_force_scale*aforce(1:3, iatom)
241 : END DO
242 : WRITE (UNIT=ounit, FMT="(T2,A,T24,3(1X,F18.12))") &
243 : "ACCINT| Total force ", tforce(1:3)
244 : END IF
245 : IF (use_virial) THEN
246 : CALL para_env%sum(avirial)
247 : IF (ounit > 0) THEN
248 : WRITE (UNIT=ounit, FMT="(/,T2,A)") "ACCINT| Debug Accurate-XCINT virial"
249 : DO i = 1, 3
250 : WRITE (UNIT=ounit, FMT="(T2,A,T24,3(1X,F18.8))") &
251 : "ACCINT| [a.u.] ", my_force_scale*avirial(i, 1:3)
252 : END DO
253 : END IF
254 : END IF
255 : END IF
256 :
257 406 : DEALLOCATE (aforce, calpha, cvalue)
258 :
259 : END IF
260 :
261 1496 : CALL timestop(handle)
262 :
263 3096 : END SUBROUTINE accint_weight_force
264 :
265 : ! **************************************************************************************************
266 : !> \brief computes the forces/virial due to atomic centered Gaussian functions
267 : !> \param e_rspace Energy density
268 : !> \param qs_env ...
269 : !> \param order ...
270 : !> \param calpha ...
271 : !> \param cvalue ...
272 : !> \param aforce ...
273 : !> \param avirial ...
274 : ! **************************************************************************************************
275 406 : SUBROUTINE gauss_grid_force(e_rspace, qs_env, order, calpha, cvalue, aforce, avirial)
276 : TYPE(pw_r3d_rs_type), INTENT(IN) :: e_rspace
277 : TYPE(qs_environment_type), POINTER :: qs_env
278 : INTEGER, INTENT(IN) :: order
279 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: calpha, cvalue
280 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: aforce
281 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: avirial
282 :
283 : CHARACTER(len=*), PARAMETER :: routineN = 'gauss_grid_force'
284 :
285 : INTEGER :: atom_a, handle, iatom, igrid, ikind, j, &
286 : natom_of_kind, ni, npme, on
287 406 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cores
288 406 : INTEGER, DIMENSION(:), POINTER :: atom_list
289 : LOGICAL :: use_virial
290 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
291 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
292 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
293 406 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
294 406 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
295 : TYPE(cell_type), POINTER :: cell
296 : TYPE(dft_control_type), POINTER :: dft_control
297 406 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
298 : TYPE(pw_env_type), POINTER :: pw_env
299 406 : TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
300 406 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_grids
301 : TYPE(realspace_grid_type), POINTER :: rs_v
302 :
303 406 : CALL timeset(routineN, handle)
304 :
305 406 : on = 2*order
306 406 : ni = ncoset(on)
307 1218 : ALLOCATE (hab(ni, 1))
308 812 : ALLOCATE (pab(ni, 1))
309 :
310 406 : NULLIFY (pw_pools, rs_grids, rs_v)
311 :
312 406 : CALL get_qs_env(qs_env, pw_env=pw_env)
313 406 : CALL pw_env_get(pw_env, pw_pools=pw_pools, rs_grids=rs_grids)
314 406 : DO igrid = 1, SIZE(pw_pools)
315 406 : IF (pw_grid_compare(e_rspace%pw_grid, pw_pools(igrid)%pool%pw_grid)) THEN
316 406 : rs_v => rs_grids(igrid)
317 406 : EXIT
318 : END IF
319 : END DO
320 406 : IF (.NOT. ASSOCIATED(rs_v)) THEN
321 0 : CPABORT("No realspace grid for Accurate-XCINT weight force")
322 : END IF
323 :
324 406 : CALL transfer_pw2rs(rs_v, e_rspace)
325 :
326 : CALL get_qs_env(qs_env, &
327 : atomic_kind_set=atomic_kind_set, &
328 : cell=cell, &
329 : dft_control=dft_control, &
330 406 : particle_set=particle_set)
331 :
332 406 : use_virial = .TRUE.
333 406 : avirial = 0.0_dp
334 5390 : aforce = 0.0_dp
335 :
336 406 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
337 :
338 1236 : DO ikind = 1, SIZE(atomic_kind_set)
339 :
340 830 : alpha = calpha(ikind)
341 830 : IF (alpha == 0.0_dp) CYCLE
342 :
343 810 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
344 :
345 810 : CALL set_polynom_coefs(pab(:, 1), order, alpha)
346 1688 : pab(:, 1) = -cvalue(ikind)*pab(:, 1)
347 :
348 2430 : ALLOCATE (cores(natom_of_kind))
349 810 : npme = 0
350 810 : cores = 0
351 :
352 2024 : DO iatom = 1, natom_of_kind
353 1214 : atom_a = atom_list(iatom)
354 1214 : ra(:) = pbc(particle_set(atom_a)%r, cell)
355 2024 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
356 : ! replicated realspace grid, split the atoms up between procs
357 1214 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
358 607 : npme = npme + 1
359 607 : cores(npme) = iatom
360 : END IF
361 : ELSE
362 0 : npme = npme + 1
363 0 : cores(npme) = iatom
364 : END IF
365 : END DO
366 :
367 1417 : DO j = 1, npme
368 :
369 607 : iatom = cores(j)
370 607 : atom_a = atom_list(iatom)
371 607 : ra(:) = pbc(particle_set(atom_a)%r, cell)
372 1282 : hab(:, 1) = 0.0_dp
373 607 : force_a(:) = 0.0_dp
374 607 : force_b(:) = 0.0_dp
375 607 : my_virial_a = 0.0_dp
376 607 : my_virial_b = 0.0_dp
377 :
378 : radius = exp_radius_very_extended(la_min=0, la_max=on, lb_min=0, lb_max=0, &
379 : ra=ra, rb=ra, rp=ra, &
380 : zetp=alpha, eps=eps_rho_rspace, &
381 : pab=pab, o1=0, o2=0, &
382 607 : prefactor=1.0_dp, cutoff=1.0_dp)
383 :
384 : CALL integrate_pgf_product(on, alpha, 0, &
385 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
386 : rs_v, hab, pab=pab, o1=0, o2=0, &
387 : radius=radius, &
388 : calculate_forces=.TRUE., force_a=force_a, &
389 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
390 607 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
391 :
392 2428 : aforce(1:3, atom_a) = aforce(1:3, atom_a) + force_a(1:3)
393 8701 : avirial = avirial + my_virial_a
394 :
395 : END DO
396 2046 : DEALLOCATE (cores)
397 :
398 : END DO
399 :
400 406 : DEALLOCATE (hab, pab)
401 :
402 406 : CALL timestop(handle)
403 :
404 812 : END SUBROUTINE gauss_grid_force
405 :
406 : ! **************************************************************************************************
407 : !> \brief computes the weight function on the PW grid
408 : !> \param rho_core ...
409 : !> \param qs_env ...
410 : !> \param order ...
411 : !> \param calpha ...
412 : !> \param ccore ...
413 : ! **************************************************************************************************
414 670 : SUBROUTINE weight_function_pwgrid(rho_core, qs_env, order, calpha, ccore)
415 : TYPE(pw_r3d_rs_type), INTENT(INOUT) :: rho_core
416 : TYPE(qs_environment_type), POINTER :: qs_env
417 : INTEGER, INTENT(IN) :: order
418 : REAL(KIND=dp), DIMENSION(:) :: calpha, ccore
419 :
420 : CHARACTER(len=*), PARAMETER :: routineN = 'weight_function_pwgrid'
421 :
422 : INTEGER :: atom_a, handle, iatom, ikind, ithread, &
423 : j, natom, ni, npme, nthread, on, &
424 : subpatch_pattern
425 670 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cores
426 670 : INTEGER, DIMENSION(:), POINTER :: atom_list
427 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
428 : REAL(KIND=dp), DIMENSION(3) :: ra
429 670 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: pab
430 670 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
431 : TYPE(cell_type), POINTER :: cell
432 : TYPE(dft_control_type), POINTER :: dft_control
433 670 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
434 : TYPE(pw_env_type), POINTER :: pw_env
435 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
436 : TYPE(pw_r3d_rs_type) :: rhoc_r
437 670 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
438 : TYPE(realspace_grid_type), POINTER :: rs_rho
439 :
440 670 : CALL timeset(routineN, handle)
441 :
442 670 : on = 2*order
443 670 : ni = ncoset(on)
444 2010 : ALLOCATE (pab(ni, 1))
445 :
446 : CALL get_qs_env(qs_env=qs_env, &
447 : atomic_kind_set=atomic_kind_set, &
448 : qs_kind_set=qs_kind_set, &
449 : cell=cell, &
450 : dft_control=dft_control, &
451 : particle_set=particle_set, &
452 670 : pw_env=pw_env)
453 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_rho, &
454 670 : auxbas_pw_pool=auxbas_pw_pool)
455 : ! be careful in parallel nsmax is chosen with multigrid in mind!
456 670 : CALL rs_grid_zero(rs_rho)
457 :
458 670 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
459 :
460 1978 : DO ikind = 1, SIZE(atomic_kind_set)
461 1308 : alpha = calpha(ikind)
462 1308 : IF (alpha == 0.0_dp) CYCLE
463 :
464 1270 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
465 :
466 1270 : CALL set_polynom_coefs(pab(:, 1), order, alpha)
467 2880 : pab(:, 1) = ccore(ikind)*pab(:, 1)
468 :
469 1270 : nthread = 1
470 1270 : ithread = 0
471 :
472 3810 : ALLOCATE (cores(natom))
473 1270 : npme = 0
474 1270 : cores = 0
475 :
476 3146 : DO iatom = 1, natom
477 3146 : IF (rs_rho%desc%parallel .AND. .NOT. rs_rho%desc%distributed) THEN
478 : ! replicated realspace grid, split the atoms up between procs
479 1876 : IF (MODULO(iatom, rs_rho%desc%group_size) == rs_rho%desc%my_pos) THEN
480 938 : npme = npme + 1
481 938 : cores(npme) = iatom
482 : END IF
483 : ELSE
484 0 : npme = npme + 1
485 0 : cores(npme) = iatom
486 : END IF
487 : END DO
488 :
489 1270 : IF (npme > 0) THEN
490 1841 : DO j = 1, npme
491 :
492 938 : iatom = cores(j)
493 938 : atom_a = atom_list(iatom)
494 938 : ra(:) = pbc(particle_set(atom_a)%r, cell)
495 938 : subpatch_pattern = 0
496 : radius = exp_radius_very_extended(la_min=0, la_max=on, &
497 : lb_min=0, lb_max=0, &
498 : ra=ra, rb=ra, rp=ra, &
499 : zetp=alpha, eps=eps_rho_rspace, &
500 : pab=pab, o1=0, o2=0, & ! without map_consistent
501 938 : prefactor=-1.0_dp, cutoff=0.0_dp)
502 :
503 : CALL collocate_pgf_product(on, alpha, 0, 0, 0.0_dp, 0, ra, &
504 : [0.0_dp, 0.0_dp, 0.0_dp], -1.0_dp, pab, 0, 0, rs_rho, &
505 : radius=radius, ga_gb_function=GRID_FUNC_AB, &
506 1841 : use_subpatch=.TRUE., subpatch_pattern=subpatch_pattern)
507 :
508 : END DO
509 : END IF
510 3248 : DEALLOCATE (cores)
511 :
512 : END DO
513 :
514 670 : DEALLOCATE (pab)
515 :
516 670 : CALL auxbas_pw_pool%create_pw(rhoc_r)
517 :
518 670 : CALL transfer_rs2pw(rs_rho, rhoc_r)
519 :
520 670 : CALL pw_transfer(rhoc_r, rho_core)
521 :
522 670 : CALL auxbas_pw_pool%give_back_pw(rhoc_r)
523 :
524 670 : CALL timestop(handle)
525 :
526 1340 : END SUBROUTINE weight_function_pwgrid
527 :
528 : ! **************************************************************************************************
529 : !> \brief Set the prefactor coefficients for the expansion of the radial weight function in
530 : !> cartesian spherical harmonics
531 : !> \param pf output coefficients
532 : !> \param order ...
533 : !> \param alpha ...
534 : ! **************************************************************************************************
535 2080 : SUBROUTINE set_polynom_coefs(pf, order, alpha)
536 : REAL(kind=dp), DIMENSION(:), INTENT(OUT) :: pf
537 : INTEGER, INTENT(IN) :: order
538 : REAL(kind=dp), INTENT(IN) :: alpha
539 :
540 : INTEGER :: i
541 : REAL(KIND=dp) :: ap
542 :
543 4568 : pf(:) = 0.0_dp
544 :
545 4184 : DO i = 0, order
546 2080 : SELECT CASE (i)
547 : CASE (0)
548 2080 : pf(1) = 1.0_dp
549 : CASE (1)
550 12 : ap = alpha
551 12 : pf(coset(2, 0, 0)) = ap
552 12 : pf(coset(0, 2, 0)) = ap
553 12 : pf(coset(0, 0, 2)) = ap
554 : CASE (2)
555 12 : ap = alpha**2/2.0_dp
556 12 : pf(coset(4, 0, 0)) = ap
557 12 : pf(coset(0, 4, 0)) = ap
558 12 : pf(coset(0, 0, 4)) = ap
559 12 : pf(coset(2, 2, 0)) = 2._dp*ap
560 12 : pf(coset(2, 0, 2)) = 2._dp*ap
561 12 : pf(coset(0, 2, 2)) = 2._dp*ap
562 : CASE (3)
563 0 : ap = alpha**3/6.0_dp
564 0 : pf(coset(6, 0, 0)) = ap
565 0 : pf(coset(0, 6, 0)) = ap
566 0 : pf(coset(0, 0, 6)) = ap
567 0 : pf(coset(4, 2, 0)) = 3._dp*ap
568 0 : pf(coset(4, 0, 2)) = 3._dp*ap
569 0 : pf(coset(2, 4, 0)) = 3._dp*ap
570 0 : pf(coset(2, 0, 4)) = 3._dp*ap
571 0 : pf(coset(0, 4, 2)) = 3._dp*ap
572 0 : pf(coset(0, 2, 4)) = 3._dp*ap
573 0 : pf(coset(2, 2, 2)) = 6._dp*ap
574 : CASE DEFAULT
575 : CALL cp_abort(__LOCATION__, &
576 : "Only 0, 1, 2, 3 are supported as the "// &
577 2104 : "polynomial order value in accuarte XC integration.")
578 : END SELECT
579 : END DO
580 :
581 2080 : END SUBROUTINE set_polynom_coefs
582 :
583 : ! **************************************************************************************************
584 : !> \brief calculates the XC density:
585 : !> order=0: exc will contain the xc energy density E_xc(r)
586 : !> order=1: exc will contain V_xc(r) * rho1(r)
587 : !> order=2: exc will contain F_xc(r) * rho1(r) * rho1(r)
588 : !> \param qs_env to get all the needed things
589 : !> \param rho_struct density
590 : !> \param rho1_struct response density
591 : !> \param order requested derivative order
592 : !> \param xc_section ...
593 : !> \param triplet ...
594 : !> \param exc Output energy density on working grid (!)
595 : !> \author JGH
596 : ! **************************************************************************************************
597 812 : SUBROUTINE xc_density(qs_env, rho_struct, rho1_struct, order, xc_section, triplet, exc)
598 :
599 : TYPE(qs_environment_type), POINTER :: qs_env
600 : TYPE(qs_rho_type), POINTER :: rho_struct, rho1_struct
601 : INTEGER, INTENT(IN) :: order
602 : TYPE(section_vals_type), POINTER :: xc_section
603 : LOGICAL, INTENT(IN) :: triplet
604 : TYPE(pw_r3d_rs_type) :: exc
605 :
606 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_density'
607 :
608 : INTEGER :: handle, ispin, myfun, nspins
609 : LOGICAL :: native_skala_grid, rho1_g_valid, &
610 : rho_g_valid, tau1_valid, tau_valid, &
611 : uf_grid
612 : REAL(KIND=dp) :: excint, factor
613 : REAL(KIND=dp), DIMENSION(3, 3) :: vdum
614 : TYPE(cell_type), POINTER :: cell
615 : TYPE(dft_control_type), POINTER :: dft_control
616 406 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
617 406 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho1_g, rho1_g_xc, rho_g, rho_g_xc
618 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
619 : TYPE(pw_env_type), POINTER :: pw_env
620 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
621 406 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho1_r_xc, rho_r, rho_r_xc, &
622 406 : tau1_r, tau1_r_xc, tau_r, tau_r_xc, &
623 406 : vxc_rho, vxc_tau
624 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights
625 406 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho0_atom_set, rho1_atom_set
626 :
627 406 : CALL timeset(routineN, handle)
628 :
629 406 : CPASSERT(ASSOCIATED(rho_struct))
630 406 : IF (order > 0) THEN
631 164 : CPASSERT(ASSOCIATED(rho1_struct))
632 : END IF
633 :
634 : ! we always get true exc (not integration weighted)
635 : CALL get_qs_env(qs_env, &
636 : dft_control=dft_control, &
637 : pw_env=pw_env, &
638 : cell=cell, &
639 : particle_set=particle_set, &
640 : rho_nlcc=rho_nlcc, &
641 406 : rho_nlcc_g=rho_nlcc_g)
642 :
643 406 : nspins = dft_control%nspins
644 406 : CPASSERT(dft_control%sic_method_id == sic_none)
645 :
646 406 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
647 406 : native_skala_grid = xc_section_uses_native_skala_grid(xc_section)
648 :
649 406 : CALL pw_zero(exc)
650 :
651 406 : IF (myfun /= xc_none) THEN
652 388 : NULLIFY (weights)
653 :
654 388 : NULLIFY (rho_r, rho_g, tau_r)
655 : CALL qs_rho_get(rho_struct, rho_r=rho_r, rho_g=rho_g, tau_r=tau_r, &
656 388 : rho_g_valid=rho_g_valid, tau_r_valid=tau_valid)
657 :
658 388 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
659 388 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
660 :
661 : ! add the nlcc densities
662 388 : IF (ASSOCIATED(rho_nlcc) .AND. order <= 1) THEN
663 8 : factor = 1.0_dp
664 16 : DO ispin = 1, nspins
665 8 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
666 16 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
667 : END DO
668 : END IF
669 :
670 388 : NULLIFY (rho_r_xc, rho_g_xc, tau_r_xc)
671 388 : IF (uf_grid .AND. order <= 1) THEN
672 36 : IF (rho_g_valid) THEN
673 36 : CALL create_density_on_pool(xc_pw_pool, rho_g, rho_r_xc, rho_g_xc)
674 : ELSE
675 0 : CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho_r, rho_r_xc, rho_g_xc)
676 : END IF
677 36 : IF (tau_valid) THEN
678 8 : CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau_r, tau_r_xc)
679 : END IF
680 : ELSE
681 352 : rho_r_xc => rho_r
682 352 : rho_g_xc => rho_g
683 352 : tau_r_xc => tau_r
684 : END IF
685 :
686 388 : NULLIFY (rho1_r_xc, rho1_g_xc, tau1_r_xc)
687 388 : IF (order >= 1) THEN
688 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, rho_g=rho1_g, tau_r=tau1_r, &
689 156 : rho_g_valid=rho1_g_valid, tau_r_valid=tau1_valid)
690 156 : IF (uf_grid) THEN
691 16 : IF (rho1_g_valid) THEN
692 16 : CALL create_density_on_pool(xc_pw_pool, rho1_g, rho1_r_xc, rho1_g_xc)
693 : ELSE
694 0 : CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, rho1_r, rho1_r_xc, rho1_g_xc)
695 : END IF
696 16 : IF (tau1_valid) THEN
697 0 : CALL create_density_on_pool_from_r(auxbas_pw_pool, xc_pw_pool, tau1_r, tau1_r_xc)
698 : END IF
699 : ELSE
700 140 : rho1_r_xc => rho1_r
701 140 : rho1_g_xc => rho1_g
702 140 : tau1_r_xc => tau1_r
703 : END IF
704 : END IF
705 :
706 388 : NULLIFY (vxc_rho, vxc_tau)
707 232 : SELECT CASE (order)
708 : CASE (0)
709 232 : IF (native_skala_grid) THEN
710 : CALL skala_gpw_weight_derivative(exc, rho_r, rho_g, tau_r, xc_section, weights, &
711 0 : xc_pw_pool, particle_set, cell)
712 : ELSE
713 232 : CALL xc_exc_pw_create(rho_r_xc, rho_g_xc, tau_r_xc, xc_section, weights, xc_pw_pool, exc)
714 : END IF
715 : CASE (1)
716 94 : IF (native_skala_grid) THEN
717 : CALL cp_abort(__LOCATION__, &
718 0 : "Native SKALA GAPW accurate-XCINT response forces are not implemented.")
719 : ELSE
720 : CALL xc_vxc_pw_create(vxc_rho, vxc_tau, exc=excint, &
721 : rho_r=rho_r_xc, rho_g=rho_g_xc, tau=tau_r_xc, &
722 : xc_section=xc_section, weights=weights, pw_pool=xc_pw_pool, &
723 94 : compute_virial=.FALSE., virial_xc=vdum)
724 : END IF
725 : CASE (2)
726 62 : IF (native_skala_grid) THEN
727 : CALL cp_abort(__LOCATION__, &
728 0 : "Native SKALA GAPW accurate-XCINT response forces are not implemented.")
729 : ELSE
730 : CALL qs_fxc_create(qs_env, rho_struct, rho1_struct, rho0_atom_set, xc_section, .FALSE., &
731 : vxc_rho, vxc_tau, rho1_atom_set, is_triplet=triplet, &
732 62 : no_weights=.TRUE., uf_grid_results=uf_grid)
733 : END IF
734 : CASE DEFAULT
735 388 : CPABORT("Derivative order not available in xc_density")
736 : END SELECT
737 :
738 388 : IF (order >= 1) THEN
739 156 : CALL pw_zero(exc)
740 156 : IF (ASSOCIATED(vxc_rho)) THEN
741 314 : DO ispin = 1, nspins
742 158 : CALL pw_multiply_with(vxc_rho(ispin), rho1_r_xc(ispin))
743 158 : CALL pw_axpy(vxc_rho(ispin), exc, 1.0_dp)
744 314 : CALL vxc_rho(ispin)%release()
745 : END DO
746 156 : DEALLOCATE (vxc_rho)
747 : END IF
748 156 : IF (ASSOCIATED(vxc_tau)) THEN
749 0 : DO ispin = 1, nspins
750 0 : CALL pw_multiply_with(vxc_tau(ispin), tau1_r_xc(ispin))
751 0 : CALL pw_axpy(vxc_tau(ispin), exc, 1.0_dp)
752 0 : CALL vxc_tau(ispin)%release()
753 : END DO
754 0 : DEALLOCATE (vxc_tau)
755 : END IF
756 : END IF
757 :
758 388 : IF (order == 2) THEN
759 62 : CALL pw_scale(exc, 0.5_dp)
760 : END IF
761 :
762 388 : IF (uf_grid .AND. order <= 1) THEN
763 36 : CALL give_back_density_on_pool(xc_pw_pool, rho_r_xc, rho_g_xc)
764 36 : IF (tau_valid) THEN
765 8 : CALL give_back_density_on_pool(xc_pw_pool, tau_r_xc)
766 : END IF
767 : END IF
768 :
769 388 : IF (uf_grid .AND. order >= 1) THEN
770 16 : CALL give_back_density_on_pool(xc_pw_pool, rho1_r_xc, rho1_g_xc)
771 16 : IF (tau1_valid) THEN
772 0 : CALL give_back_density_on_pool(xc_pw_pool, tau1_r_xc)
773 : END IF
774 : END IF
775 :
776 : ! romove the nlcc densities
777 388 : IF (ASSOCIATED(rho_nlcc) .AND. order <= 1) THEN
778 8 : factor = -1.0_dp
779 16 : DO ispin = 1, nspins
780 8 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
781 16 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
782 : END DO
783 : END IF
784 :
785 : END IF
786 :
787 406 : CALL timestop(handle)
788 :
789 406 : END SUBROUTINE xc_density
790 : ! **************************************************************************************************
791 : !> \brief transfers a g-space density to a given PW pool and creates its r-space representation
792 : !> \param pw_pool ...
793 : !> \param rho_g_in ...
794 : !> \param rho_r_out ...
795 : !> \param rho_g_out ...
796 : ! **************************************************************************************************
797 52 : SUBROUTINE create_density_on_pool(pw_pool, rho_g_in, rho_r_out, rho_g_out)
798 : TYPE(pw_pool_type), POINTER :: pw_pool
799 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_in
800 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_out
801 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g_out
802 :
803 : INTEGER :: ispin, nspins
804 :
805 52 : CPASSERT(ASSOCIATED(pw_pool))
806 52 : CPASSERT(ASSOCIATED(rho_g_in))
807 :
808 52 : nspins = SIZE(rho_g_in)
809 364 : ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
810 104 : DO ispin = 1, nspins
811 52 : CALL pw_pool%create_pw(rho_g_out(ispin))
812 52 : CALL pw_pool%create_pw(rho_r_out(ispin))
813 52 : CALL pw_transfer(rho_g_in(ispin), rho_g_out(ispin))
814 104 : CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
815 : END DO
816 :
817 52 : END SUBROUTINE create_density_on_pool
818 :
819 : ! **************************************************************************************************
820 : !> \brief transfers an r-space density to a given PW pool and creates its g-space representation
821 : !> \param source_pw_pool ...
822 : !> \param target_pw_pool ...
823 : !> \param rho_r_in ...
824 : !> \param rho_r_out ...
825 : !> \param rho_g_out ...
826 : ! **************************************************************************************************
827 8 : SUBROUTINE create_density_on_pool_from_r(source_pw_pool, target_pw_pool, rho_r_in, rho_r_out, rho_g_out)
828 : TYPE(pw_pool_type), POINTER :: source_pw_pool, target_pw_pool
829 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r_in, rho_r_out
830 : TYPE(pw_c1d_gs_type), DIMENSION(:), OPTIONAL, &
831 : POINTER :: rho_g_out
832 :
833 : INTEGER :: ispin, nspins
834 : TYPE(pw_c1d_gs_type) :: rho_g_aux, rho_g_in
835 :
836 0 : CPASSERT(ASSOCIATED(source_pw_pool))
837 8 : CPASSERT(ASSOCIATED(target_pw_pool))
838 8 : CPASSERT(ASSOCIATED(rho_r_in))
839 :
840 8 : nspins = SIZE(rho_r_in)
841 8 : IF (PRESENT(rho_g_out)) THEN
842 0 : ALLOCATE (rho_r_out(nspins), rho_g_out(nspins))
843 0 : DO ispin = 1, nspins
844 0 : CALL source_pw_pool%create_pw(rho_g_in)
845 0 : CALL target_pw_pool%create_pw(rho_g_out(ispin))
846 0 : CALL target_pw_pool%create_pw(rho_r_out(ispin))
847 0 : CALL pw_transfer(rho_r_in(ispin), rho_g_in)
848 0 : CALL pw_transfer(rho_g_in, rho_g_out(ispin))
849 0 : CALL pw_transfer(rho_g_out(ispin), rho_r_out(ispin))
850 0 : CALL source_pw_pool%give_back_pw(rho_g_in)
851 : END DO
852 : ELSE
853 32 : ALLOCATE (rho_r_out(nspins))
854 16 : DO ispin = 1, nspins
855 8 : CALL source_pw_pool%create_pw(rho_g_in)
856 8 : CALL target_pw_pool%create_pw(rho_g_aux)
857 8 : CALL target_pw_pool%create_pw(rho_r_out(ispin))
858 8 : CALL pw_transfer(rho_r_in(ispin), rho_g_in)
859 8 : CALL pw_transfer(rho_g_in, rho_g_aux)
860 8 : CALL pw_transfer(rho_g_aux, rho_r_out(ispin))
861 8 : CALL source_pw_pool%give_back_pw(rho_g_in)
862 16 : CALL source_pw_pool%give_back_pw(rho_g_aux)
863 : END DO
864 : END IF
865 :
866 8 : END SUBROUTINE create_density_on_pool_from_r
867 :
868 : ! **************************************************************************************************
869 : !> \brief returns temporary density arrays to the given PW pool
870 : !> \param pw_pool ...
871 : !> \param rho_r ...
872 : !> \param rho_g ...
873 : ! **************************************************************************************************
874 60 : SUBROUTINE give_back_density_on_pool(pw_pool, rho_r, rho_g)
875 : TYPE(pw_pool_type), POINTER :: pw_pool
876 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
877 : TYPE(pw_c1d_gs_type), DIMENSION(:), OPTIONAL, &
878 : POINTER :: rho_g
879 :
880 : INTEGER :: ispin
881 :
882 60 : CPASSERT(ASSOCIATED(pw_pool))
883 :
884 60 : IF (ASSOCIATED(rho_r)) THEN
885 120 : DO ispin = 1, SIZE(rho_r)
886 120 : CALL pw_pool%give_back_pw(rho_r(ispin))
887 : END DO
888 60 : DEALLOCATE (rho_r)
889 : END IF
890 60 : IF (PRESENT(rho_g)) THEN
891 52 : IF (ASSOCIATED(rho_g)) THEN
892 104 : DO ispin = 1, SIZE(rho_g)
893 104 : CALL pw_pool%give_back_pw(rho_g(ispin))
894 : END DO
895 52 : DEALLOCATE (rho_g)
896 : END IF
897 : END IF
898 :
899 60 : END SUBROUTINE give_back_density_on_pool
900 :
901 : END MODULE accint_weights_forces
|
