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