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_zero
34 : USE pw_pool_types, ONLY: pw_pool_type
35 : USE pw_types, ONLY: pw_c1d_gs_type,&
36 : pw_r3d_rs_type
37 : USE qs_environment_types, ONLY: get_qs_env,&
38 : qs_environment_type
39 : USE qs_force_types, ONLY: qs_force_type
40 : USE qs_fxc, ONLY: qs_fxc_analytic
41 : USE qs_ks_types, ONLY: get_ks_env,&
42 : qs_ks_env_type
43 : USE qs_rho_types, ONLY: qs_rho_get,&
44 : qs_rho_type
45 : USE realspace_grid_types, ONLY: realspace_grid_type,&
46 : transfer_pw2rs
47 : USE virial_types, ONLY: virial_type
48 : USE xc, ONLY: xc_exc_pw_create,&
49 : xc_vxc_pw_create
50 : #include "./base/base_uses.f90"
51 :
52 : IMPLICIT NONE
53 :
54 : PRIVATE
55 :
56 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
57 :
58 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'accint_weights_forces'
59 :
60 : PUBLIC :: accint_weight_force
61 :
62 : CONTAINS
63 :
64 : ! **************************************************************************************************
65 : !> \brief ...
66 : !> \param qs_env ...
67 : !> \param rho ...
68 : !> \param rho1 ...
69 : !> \param order ...
70 : !> \param xc_section ...
71 : !> \param triplet ...
72 : ! **************************************************************************************************
73 6577 : SUBROUTINE accint_weight_force(qs_env, rho, rho1, order, xc_section, triplet)
74 : TYPE(qs_environment_type), POINTER :: qs_env
75 : TYPE(qs_rho_type), POINTER :: rho, rho1
76 : INTEGER, INTENT(IN) :: order
77 : TYPE(section_vals_type), POINTER :: xc_section
78 : LOGICAL, INTENT(IN), OPTIONAL :: triplet
79 :
80 : CHARACTER(len=*), PARAMETER :: routineN = 'accint_weight_force'
81 :
82 : INTEGER :: atom_a, handle, iatom, ikind, natom, &
83 : natom_of_kind, nkind
84 6577 : INTEGER, DIMENSION(:), POINTER :: atom_list
85 : LOGICAL :: lr_triplet, use_virial
86 6577 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: calpha, cvalue
87 6577 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: aforce
88 : REAL(KIND=dp), DIMENSION(3, 3) :: avirial
89 6577 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
90 : TYPE(dft_control_type), POINTER :: dft_control
91 : TYPE(pw_env_type), POINTER :: pw_env
92 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
93 : TYPE(pw_r3d_rs_type) :: e_rspace
94 6577 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
95 : TYPE(qs_ks_env_type), POINTER :: ks_env
96 : TYPE(virial_type), POINTER :: virial
97 :
98 6577 : CALL timeset(routineN, handle)
99 :
100 6577 : CALL get_qs_env(qs_env, dft_control=dft_control)
101 :
102 6577 : IF (dft_control%qs_control%gapw_control%accurate_xcint) THEN
103 :
104 278 : CALL get_qs_env(qs_env=qs_env, force=force, virial=virial)
105 278 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
106 :
107 278 : CALL get_qs_env(qs_env, natom=natom, nkind=nkind)
108 834 : ALLOCATE (aforce(3, natom))
109 1390 : ALLOCATE (calpha(nkind), cvalue(nkind))
110 896 : cvalue = 1.0_dp
111 896 : calpha(1:nkind) = dft_control%qs_control%gapw_control%aw(1:nkind)
112 :
113 278 : CALL get_qs_env(qs_env, ks_env=ks_env, pw_env=pw_env)
114 278 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
115 278 : CALL auxbas_pw_pool%create_pw(e_rspace)
116 :
117 278 : lr_triplet = .FALSE.
118 278 : IF (PRESENT(triplet)) lr_triplet = triplet
119 :
120 278 : CALL xc_density(ks_env, rho, rho1, order, xc_section, lr_triplet, e_rspace)
121 278 : CALL pw_scale(e_rspace, e_rspace%pw_grid%dvol)
122 :
123 278 : CALL gauss_grid_force(e_rspace, qs_env, calpha, cvalue, aforce, avirial)
124 :
125 278 : CALL auxbas_pw_pool%give_back_pw(e_rspace)
126 :
127 278 : CALL get_qs_env(qs_env, atomic_kind_set=atomic_kind_set)
128 896 : DO ikind = 1, nkind
129 618 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
130 1734 : DO iatom = 1, natom_of_kind
131 838 : atom_a = atom_list(iatom)
132 : force(ikind)%rho_elec(1:3, iatom) = &
133 3970 : force(ikind)%rho_elec(1:3, iatom) + aforce(1:3, atom_a)
134 : END DO
135 : END DO
136 278 : IF (use_virial) THEN
137 0 : virial%pv_exc = virial%pv_exc + avirial
138 0 : virial%pv_virial = virial%pv_virial + avirial
139 : END IF
140 :
141 278 : DEALLOCATE (aforce)
142 :
143 : END IF
144 :
145 6577 : CALL timestop(handle)
146 :
147 13154 : END SUBROUTINE accint_weight_force
148 :
149 : ! **************************************************************************************************
150 : !> \brief computes the forces/virial due to atomic centered Gaussian functions
151 : !> \param e_rspace Energy density
152 : !> \param qs_env ...
153 : !> \param calpha ...
154 : !> \param cvalue ...
155 : !> \param aforce ...
156 : !> \param avirial ...
157 : ! **************************************************************************************************
158 278 : SUBROUTINE gauss_grid_force(e_rspace, qs_env, calpha, cvalue, aforce, avirial)
159 : TYPE(pw_r3d_rs_type), INTENT(IN) :: e_rspace
160 : TYPE(qs_environment_type), POINTER :: qs_env
161 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: calpha, cvalue
162 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: aforce
163 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT) :: avirial
164 :
165 : CHARACTER(len=*), PARAMETER :: routineN = 'gauss_grid_force'
166 :
167 : INTEGER :: atom_a, handle, iatom, ikind, j, &
168 : natom_of_kind, npme
169 278 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
170 : LOGICAL :: use_virial
171 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
172 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
173 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
174 278 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
175 278 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
176 : TYPE(cell_type), POINTER :: cell
177 : TYPE(dft_control_type), POINTER :: dft_control
178 278 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
179 : TYPE(pw_env_type), POINTER :: pw_env
180 : TYPE(realspace_grid_type), POINTER :: rs_v
181 :
182 278 : CALL timeset(routineN, handle)
183 :
184 278 : ALLOCATE (cores(1))
185 278 : ALLOCATE (hab(1, 1))
186 278 : ALLOCATE (pab(1, 1))
187 :
188 278 : CALL get_qs_env(qs_env, pw_env=pw_env)
189 278 : CALL pw_env_get(pw_env, auxbas_rs_grid=rs_v)
190 :
191 278 : CALL transfer_pw2rs(rs_v, e_rspace)
192 :
193 : CALL get_qs_env(qs_env, &
194 : atomic_kind_set=atomic_kind_set, &
195 : cell=cell, &
196 : dft_control=dft_control, &
197 278 : particle_set=particle_set)
198 :
199 278 : use_virial = .TRUE.
200 278 : avirial = 0.0_dp
201 3630 : aforce = 0.0_dp
202 :
203 278 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
204 :
205 896 : DO ikind = 1, SIZE(atomic_kind_set)
206 :
207 618 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
208 :
209 618 : alpha = calpha(ikind)
210 618 : pab(1, 1) = -cvalue(ikind)
211 618 : IF (alpha == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
212 :
213 606 : CALL reallocate(cores, 1, natom_of_kind)
214 606 : npme = 0
215 1432 : cores = 0
216 :
217 1432 : DO iatom = 1, natom_of_kind
218 826 : atom_a = atom_list(iatom)
219 826 : ra(:) = pbc(particle_set(atom_a)%r, cell)
220 1432 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
221 : ! replicated realspace grid, split the atoms up between procs
222 826 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
223 413 : npme = npme + 1
224 413 : cores(npme) = iatom
225 : END IF
226 : ELSE
227 0 : npme = npme + 1
228 0 : cores(npme) = iatom
229 : END IF
230 : END DO
231 :
232 1915 : DO j = 1, npme
233 :
234 413 : iatom = cores(j)
235 413 : atom_a = atom_list(iatom)
236 413 : ra(:) = pbc(particle_set(atom_a)%r, cell)
237 413 : hab(1, 1) = 0.0_dp
238 413 : force_a(:) = 0.0_dp
239 413 : force_b(:) = 0.0_dp
240 413 : my_virial_a = 0.0_dp
241 413 : my_virial_b = 0.0_dp
242 :
243 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
244 : ra=ra, rb=ra, rp=ra, &
245 : zetp=alpha, eps=eps_rho_rspace, &
246 : pab=pab, o1=0, o2=0, &
247 413 : prefactor=1.0_dp, cutoff=1.0_dp)
248 :
249 : CALL integrate_pgf_product(0, alpha, 0, &
250 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
251 : rs_v, hab, pab=pab, o1=0, o2=0, &
252 : radius=radius, &
253 : calculate_forces=.TRUE., force_a=force_a, &
254 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
255 413 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
256 :
257 1652 : aforce(1:3, atom_a) = aforce(1:3, atom_a) + force_a(1:3)
258 5987 : avirial = avirial + my_virial_a
259 :
260 : END DO
261 :
262 : END DO
263 :
264 278 : DEALLOCATE (hab, pab, cores)
265 :
266 278 : CALL timestop(handle)
267 :
268 278 : END SUBROUTINE gauss_grid_force
269 :
270 : ! **************************************************************************************************
271 : !> \brief calculates the XC density:
272 : !> order=0: exc will contain the xc energy density E_xc(r)
273 : !> order=1: exc will contain V_xc(r) * rho1(r)
274 : !> order=2: exc will contain F_xc(r) * rho1(r) * rho1(r)
275 : !> \param ks_env to get all the needed things
276 : !> \param rho_struct density
277 : !> \param rho1_struct response density
278 : !> \param order requested derivative order
279 : !> \param xc_section ...
280 : !> \param triplet ...
281 : !> \param exc ...
282 : !> \author JGH
283 : ! **************************************************************************************************
284 556 : SUBROUTINE xc_density(ks_env, rho_struct, rho1_struct, order, xc_section, triplet, exc)
285 :
286 : TYPE(qs_ks_env_type), POINTER :: ks_env
287 : TYPE(qs_rho_type), POINTER :: rho_struct, rho1_struct
288 : INTEGER, INTENT(IN) :: order
289 : TYPE(section_vals_type), POINTER :: xc_section
290 : LOGICAL, INTENT(IN) :: triplet
291 : TYPE(pw_r3d_rs_type) :: exc
292 :
293 : CHARACTER(len=*), PARAMETER :: routineN = 'xc_density'
294 :
295 : INTEGER :: handle, ispin, myfun, nspins
296 : LOGICAL :: uf_grid
297 : REAL(KIND=dp) :: excint, factor
298 : REAL(KIND=dp), DIMENSION(3, 3) :: vdum
299 : TYPE(cell_type), POINTER :: cell
300 : TYPE(dft_control_type), POINTER :: dft_control
301 278 : TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
302 : TYPE(pw_c1d_gs_type), POINTER :: rho_nlcc_g
303 : TYPE(pw_env_type), POINTER :: pw_env
304 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool, xc_pw_pool
305 278 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho_r, tau1_r, tau_r, vxc_rho, &
306 278 : vxc_tau
307 : TYPE(pw_r3d_rs_type), POINTER :: rho_nlcc, weights
308 :
309 278 : CALL timeset(routineN, handle)
310 :
311 : ! we always get true exc (not integration weighted)
312 278 : NULLIFY (weights)
313 :
314 : CALL get_ks_env(ks_env, &
315 : dft_control=dft_control, &
316 : pw_env=pw_env, &
317 : cell=cell, &
318 : rho_nlcc=rho_nlcc, &
319 278 : rho_nlcc_g=rho_nlcc_g)
320 :
321 278 : CALL qs_rho_get(rho_struct, rho_r=rho_r, rho_g=rho_g, tau_r=tau_r)
322 :
323 278 : nspins = dft_control%nspins
324 :
325 278 : CALL section_vals_val_get(xc_section, "XC_FUNCTIONAL%_SECTION_PARAMETERS_", i_val=myfun)
326 :
327 278 : CALL pw_env_get(pw_env, xc_pw_pool=xc_pw_pool, auxbas_pw_pool=auxbas_pw_pool)
328 278 : uf_grid = .NOT. pw_grid_compare(auxbas_pw_pool%pw_grid, xc_pw_pool%pw_grid)
329 278 : IF (uf_grid) THEN
330 0 : CALL cp_warn(__LOCATION__, "Fine grid option not possible with energy density")
331 0 : CPABORT("Fine Grid in xc_density")
332 : END IF
333 :
334 278 : CALL pw_zero(exc)
335 :
336 278 : IF (myfun /= xc_none) THEN
337 :
338 266 : CPASSERT(ASSOCIATED(rho_struct))
339 266 : CPASSERT(dft_control%sic_method_id == sic_none)
340 :
341 : ! add the nlcc densities
342 266 : IF (ASSOCIATED(rho_nlcc) .AND. order <= 1) THEN
343 6 : factor = 1.0_dp
344 12 : DO ispin = 1, nspins
345 6 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
346 12 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
347 : END DO
348 : END IF
349 :
350 266 : NULLIFY (vxc_rho, vxc_tau)
351 394 : SELECT CASE (order)
352 : CASE (0)
353 : ! we could reduce to energy only here
354 128 : CALL xc_exc_pw_create(rho_r, rho_g, tau_r, xc_section, weights, xc_pw_pool, exc)
355 : CASE (1)
356 84 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, tau_r=tau1_r)
357 : CALL xc_vxc_pw_create(vxc_rho=vxc_rho, vxc_tau=vxc_tau, rho_r=rho_r, &
358 : rho_g=rho_g, tau=tau_r, exc=excint, &
359 : xc_section=xc_section, &
360 : weights=weights, pw_pool=xc_pw_pool, &
361 : compute_virial=.FALSE., &
362 84 : virial_xc=vdum)
363 : CASE (2)
364 54 : CALL qs_rho_get(rho1_struct, rho_r=rho1_r, tau_r=tau1_r)
365 : CALL qs_fxc_analytic(rho_struct, rho1_r, tau1_r, xc_section, weights, xc_pw_pool, &
366 54 : triplet, vxc_rho, vxc_tau)
367 : CASE DEFAULT
368 266 : CPABORT("Derivative order not available in xc_density")
369 : END SELECT
370 :
371 : ! remove the nlcc densities (keep stuff in original state)
372 266 : IF (ASSOCIATED(rho_nlcc) .AND. order <= 1) THEN
373 6 : factor = -1.0_dp
374 12 : DO ispin = 1, dft_control%nspins
375 6 : CALL pw_axpy(rho_nlcc, rho_r(ispin), factor)
376 12 : CALL pw_axpy(rho_nlcc_g, rho_g(ispin), factor)
377 : END DO
378 : END IF
379 : !
380 138 : SELECT CASE (order)
381 : CASE (0)
382 : !
383 : CASE (1, 2)
384 138 : CALL pw_zero(exc)
385 138 : IF (ASSOCIATED(vxc_rho)) THEN
386 278 : DO ispin = 1, nspins
387 140 : CALL pw_multiply_with(vxc_rho(ispin), rho1_r(ispin))
388 140 : CALL pw_axpy(vxc_rho(ispin), exc, 1.0_dp)
389 278 : CALL vxc_rho(ispin)%release()
390 : END DO
391 138 : DEALLOCATE (vxc_rho)
392 : END IF
393 138 : IF (ASSOCIATED(vxc_tau)) THEN
394 0 : DO ispin = 1, nspins
395 0 : CALL pw_multiply_with(vxc_tau(ispin), tau1_r(ispin))
396 0 : CALL pw_axpy(vxc_tau(ispin), exc, 1.0_dp)
397 0 : CALL vxc_tau(ispin)%release()
398 : END DO
399 0 : DEALLOCATE (vxc_tau)
400 : END IF
401 : CASE DEFAULT
402 266 : CPABORT("Derivative order not available in xc_density")
403 : END SELECT
404 :
405 266 : IF (order == 2) THEN
406 54 : CALL pw_scale(exc, 0.5_dp)
407 : END IF
408 :
409 : END IF
410 :
411 278 : CALL timestop(handle)
412 :
413 278 : END SUBROUTINE xc_density
414 :
415 : END MODULE accint_weights_forces
|
