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