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