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