Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Build up the plane wave density by collocating the primitive Gaussian
10 : !> functions (pgf).
11 : !> \par History
12 : !> Joost VandeVondele (02.2002)
13 : !> 1) rewrote collocate_pgf for increased accuracy and speed
14 : !> 2) collocate_core hack for PGI compiler
15 : !> 3) added multiple grid feature
16 : !> 4) new way to go over the grid
17 : !> Joost VandeVondele (05.2002)
18 : !> 1) prelim. introduction of the real space grid type
19 : !> JGH [30.08.02] multigrid arrays independent from potential
20 : !> JGH [17.07.03] distributed real space code
21 : !> JGH [23.11.03] refactoring and new loop ordering
22 : !> JGH [04.12.03] OpneMP parallelization of main loops
23 : !> Joost VandeVondele (12.2003)
24 : !> 1) modified to compute tau
25 : !> Joost removed incremental build feature
26 : !> Joost introduced map consistent
27 : !> Rewrote grid integration/collocation routines, [Joost VandeVondele,03.2007]
28 : !> \author Matthias Krack (03.04.2001)
29 : ! **************************************************************************************************
30 : MODULE qs_integrate_potential_single
31 : USE ao_util, ONLY: exp_radius_very_extended
32 : USE atomic_kind_types, ONLY: atomic_kind_type,&
33 : get_atomic_kind,&
34 : get_atomic_kind_set
35 : USE atprop_types, ONLY: atprop_array_init,&
36 : atprop_type
37 : USE basis_set_types, ONLY: get_gto_basis_set,&
38 : gto_basis_set_type
39 : USE cell_types, ONLY: cell_type,&
40 : pbc
41 : USE cp_control_types, ONLY: dft_control_type
42 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
43 : dbcsr_type
44 : USE external_potential_types, ONLY: get_potential,&
45 : gth_potential_type
46 : USE gaussian_gridlevels, ONLY: gaussian_gridlevel,&
47 : gridlevel_info_type
48 : USE grid_api, ONLY: integrate_pgf_product
49 : USE kinds, ONLY: dp
50 : USE lri_environment_types, ONLY: lri_kind_type
51 : USE memory_utilities, ONLY: reallocate
52 : USE message_passing, ONLY: mp_para_env_type
53 : USE orbital_pointers, ONLY: coset,&
54 : ncoset
55 : USE particle_types, ONLY: particle_type
56 : USE pw_env_types, ONLY: pw_env_get,&
57 : pw_env_type
58 : USE pw_types, ONLY: pw_r3d_rs_type
59 : USE qs_environment_types, ONLY: get_qs_env,&
60 : qs_environment_type
61 : USE qs_force_types, ONLY: qs_force_type
62 : USE qs_kind_types, ONLY: get_qs_kind,&
63 : get_qs_kind_set,&
64 : qs_kind_type
65 : USE realspace_grid_types, ONLY: map_gaussian_here,&
66 : realspace_grid_type,&
67 : rs_grid_zero,&
68 : transfer_pw2rs
69 : USE rs_pw_interface, ONLY: potential_pw2rs
70 : USE virial_types, ONLY: virial_type
71 : #include "./base/base_uses.f90"
72 :
73 : IMPLICIT NONE
74 :
75 : PRIVATE
76 :
77 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
78 :
79 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_integrate_potential_single'
80 :
81 : ! *** Public subroutines ***
82 : ! *** Don't include this routines directly, use the interface to
83 : ! *** qs_integrate_potential
84 :
85 : PUBLIC :: integrate_v_rspace_one_center, &
86 : integrate_v_rspace_diagonal, &
87 : integrate_v_core_rspace, &
88 : integrate_v_gaussian_rspace, &
89 : integrate_function, &
90 : integrate_ppl_rspace, &
91 : integrate_rho_nlcc
92 :
93 : CONTAINS
94 :
95 : ! **************************************************************************************************
96 : !> \brief computes the forces/virial due to the local pseudopotential
97 : !> \param rho_rspace ...
98 : !> \param qs_env ...
99 : ! **************************************************************************************************
100 14 : SUBROUTINE integrate_ppl_rspace(rho_rspace, qs_env)
101 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_rspace
102 : TYPE(qs_environment_type), POINTER :: qs_env
103 :
104 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_ppl_rspace'
105 :
106 : INTEGER :: atom_a, handle, iatom, ikind, j, lppl, &
107 : n, natom_of_kind, ni, npme
108 14 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
109 : LOGICAL :: use_virial
110 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
111 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
112 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
113 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl
114 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
115 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
116 : TYPE(cell_type), POINTER :: cell
117 : TYPE(dft_control_type), POINTER :: dft_control
118 : TYPE(gth_potential_type), POINTER :: gth_potential
119 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
120 : TYPE(pw_env_type), POINTER :: pw_env
121 14 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
122 14 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
123 : TYPE(realspace_grid_type), POINTER :: rs_v
124 : TYPE(virial_type), POINTER :: virial
125 :
126 14 : CALL timeset(routineN, handle)
127 :
128 14 : NULLIFY (pw_env, cores)
129 :
130 14 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
131 14 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
132 :
133 14 : CALL transfer_pw2rs(rs_v, rho_rspace)
134 :
135 : CALL get_qs_env(qs_env=qs_env, &
136 : atomic_kind_set=atomic_kind_set, &
137 : qs_kind_set=qs_kind_set, &
138 : cell=cell, &
139 : dft_control=dft_control, &
140 : particle_set=particle_set, &
141 : pw_env=pw_env, &
142 14 : force=force, virial=virial)
143 :
144 14 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
145 :
146 14 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
147 :
148 34 : DO ikind = 1, SIZE(atomic_kind_set)
149 :
150 20 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
151 20 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
152 :
153 20 : IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
154 20 : CALL get_potential(potential=gth_potential, alpha_ppl=alpha, nexp_ppl=lppl, cexp_ppl=cexp_ppl)
155 :
156 20 : IF (lppl <= 0) CYCLE
157 :
158 20 : ni = ncoset(2*lppl - 2)
159 100 : ALLOCATE (hab(ni, 1), pab(ni, 1))
160 240 : pab = 0._dp
161 :
162 20 : CALL reallocate(cores, 1, natom_of_kind)
163 20 : npme = 0
164 70 : cores = 0
165 :
166 : ! prepare core function
167 60 : DO j = 1, lppl
168 20 : SELECT CASE (j)
169 : CASE (1)
170 20 : pab(1, 1) = cexp_ppl(1)
171 : CASE (2)
172 20 : n = coset(2, 0, 0)
173 20 : pab(n, 1) = cexp_ppl(2)
174 20 : n = coset(0, 2, 0)
175 20 : pab(n, 1) = cexp_ppl(2)
176 20 : n = coset(0, 0, 2)
177 20 : pab(n, 1) = cexp_ppl(2)
178 : CASE (3)
179 0 : n = coset(4, 0, 0)
180 0 : pab(n, 1) = cexp_ppl(3)
181 0 : n = coset(0, 4, 0)
182 0 : pab(n, 1) = cexp_ppl(3)
183 0 : n = coset(0, 0, 4)
184 0 : pab(n, 1) = cexp_ppl(3)
185 0 : n = coset(2, 2, 0)
186 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
187 0 : n = coset(2, 0, 2)
188 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
189 0 : n = coset(0, 2, 2)
190 0 : pab(n, 1) = 2._dp*cexp_ppl(3)
191 : CASE (4)
192 0 : n = coset(6, 0, 0)
193 0 : pab(n, 1) = cexp_ppl(4)
194 0 : n = coset(0, 6, 0)
195 0 : pab(n, 1) = cexp_ppl(4)
196 0 : n = coset(0, 0, 6)
197 0 : pab(n, 1) = cexp_ppl(4)
198 0 : n = coset(4, 2, 0)
199 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
200 0 : n = coset(4, 0, 2)
201 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
202 0 : n = coset(2, 4, 0)
203 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
204 0 : n = coset(2, 0, 4)
205 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
206 0 : n = coset(0, 4, 2)
207 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
208 0 : n = coset(0, 2, 4)
209 0 : pab(n, 1) = 3._dp*cexp_ppl(4)
210 0 : n = coset(2, 2, 2)
211 0 : pab(n, 1) = 6._dp*cexp_ppl(4)
212 : CASE DEFAULT
213 40 : CPABORT("")
214 : END SELECT
215 : END DO
216 :
217 70 : DO iatom = 1, natom_of_kind
218 50 : atom_a = atom_list(iatom)
219 50 : ra(:) = pbc(particle_set(atom_a)%r, cell)
220 70 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
221 : ! replicated realspace grid, split the atoms up between procs
222 50 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
223 25 : npme = npme + 1
224 25 : 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 45 : DO j = 1, npme
233 :
234 25 : iatom = cores(j)
235 25 : atom_a = atom_list(iatom)
236 25 : ra(:) = pbc(particle_set(atom_a)%r, cell)
237 275 : hab(:, 1) = 0.0_dp
238 25 : force_a(:) = 0.0_dp
239 25 : force_b(:) = 0.0_dp
240 25 : IF (use_virial) THEN
241 0 : my_virial_a = 0.0_dp
242 0 : my_virial_b = 0.0_dp
243 : END IF
244 25 : ni = 2*lppl - 2
245 :
246 : radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
247 : ra=ra, rb=ra, rp=ra, &
248 : zetp=alpha, eps=eps_rho_rspace, &
249 : pab=pab, o1=0, o2=0, & ! without map_consistent
250 25 : prefactor=1.0_dp, cutoff=1.0_dp)
251 :
252 : CALL integrate_pgf_product(ni, alpha, 0, &
253 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
254 : rs_v, hab, pab=pab, o1=0, o2=0, &
255 : radius=radius, &
256 : calculate_forces=.TRUE., force_a=force_a, &
257 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
258 25 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
259 :
260 : force(ikind)%gth_ppl(:, iatom) = &
261 100 : force(ikind)%gth_ppl(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
262 :
263 45 : IF (use_virial) THEN
264 0 : virial%pv_ppl = virial%pv_ppl + my_virial_a*rho_rspace%pw_grid%dvol
265 0 : virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
266 0 : CPABORT("Virial not debuged for CORE_PPL")
267 : END IF
268 : END DO
269 :
270 74 : DEALLOCATE (hab, pab)
271 :
272 : END DO
273 :
274 14 : DEALLOCATE (cores)
275 :
276 14 : CALL timestop(handle)
277 :
278 14 : END SUBROUTINE integrate_ppl_rspace
279 :
280 : ! **************************************************************************************************
281 : !> \brief computes the forces/virial due to the nlcc pseudopotential
282 : !> \param rho_rspace ...
283 : !> \param qs_env ...
284 : ! **************************************************************************************************
285 30 : SUBROUTINE integrate_rho_nlcc(rho_rspace, qs_env)
286 : TYPE(pw_r3d_rs_type), INTENT(IN) :: rho_rspace
287 : TYPE(qs_environment_type), POINTER :: qs_env
288 :
289 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_rho_nlcc'
290 :
291 : INTEGER :: atom_a, handle, iatom, iexp_nlcc, ikind, &
292 : ithread, j, n, natom, nc, nexp_nlcc, &
293 : ni, npme, nthread
294 30 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores, nct_nlcc
295 : LOGICAL :: nlcc, use_virial
296 : REAL(KIND=dp) :: alpha, eps_rho_rspace, radius
297 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
298 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
299 30 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha_nlcc
300 30 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cval_nlcc, hab, pab
301 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
302 : TYPE(cell_type), POINTER :: cell
303 : TYPE(dft_control_type), POINTER :: dft_control
304 : TYPE(gth_potential_type), POINTER :: gth_potential
305 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
306 : TYPE(pw_env_type), POINTER :: pw_env
307 30 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
308 30 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
309 : TYPE(realspace_grid_type), POINTER :: rs_v
310 : TYPE(virial_type), POINTER :: virial
311 :
312 30 : CALL timeset(routineN, handle)
313 :
314 30 : NULLIFY (pw_env, cores)
315 :
316 30 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
317 30 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
318 :
319 30 : CALL transfer_pw2rs(rs_v, rho_rspace)
320 :
321 : CALL get_qs_env(qs_env=qs_env, &
322 : atomic_kind_set=atomic_kind_set, &
323 : qs_kind_set=qs_kind_set, &
324 : cell=cell, &
325 : dft_control=dft_control, &
326 : particle_set=particle_set, &
327 : pw_env=pw_env, &
328 30 : force=force, virial=virial)
329 :
330 30 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
331 :
332 30 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
333 :
334 74 : DO ikind = 1, SIZE(atomic_kind_set)
335 :
336 44 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
337 44 : CALL get_qs_kind(qs_kind_set(ikind), gth_potential=gth_potential)
338 :
339 44 : IF (.NOT. ASSOCIATED(gth_potential)) CYCLE
340 : CALL get_potential(potential=gth_potential, nlcc_present=nlcc, nexp_nlcc=nexp_nlcc, &
341 44 : alpha_nlcc=alpha_nlcc, nct_nlcc=nct_nlcc, cval_nlcc=cval_nlcc)
342 :
343 44 : IF (.NOT. nlcc) CYCLE
344 :
345 202 : DO iexp_nlcc = 1, nexp_nlcc
346 :
347 42 : alpha = alpha_nlcc(iexp_nlcc)
348 42 : nc = nct_nlcc(iexp_nlcc)
349 :
350 42 : ni = ncoset(2*nc - 2)
351 :
352 42 : nthread = 1
353 42 : ithread = 0
354 :
355 168 : ALLOCATE (hab(ni, 1), pab(ni, 1))
356 270 : pab = 0._dp
357 :
358 42 : CALL reallocate(cores, 1, natom)
359 42 : npme = 0
360 172 : cores = 0
361 :
362 : ! prepare core function
363 100 : DO j = 1, nc
364 42 : SELECT CASE (j)
365 : CASE (1)
366 42 : pab(1, 1) = cval_nlcc(1, iexp_nlcc)
367 : CASE (2)
368 16 : n = coset(2, 0, 0)
369 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
370 16 : n = coset(0, 2, 0)
371 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
372 16 : n = coset(0, 0, 2)
373 16 : pab(n, 1) = cval_nlcc(2, iexp_nlcc)/alpha**2
374 : CASE (3)
375 0 : n = coset(4, 0, 0)
376 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
377 0 : n = coset(0, 4, 0)
378 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
379 0 : n = coset(0, 0, 4)
380 0 : pab(n, 1) = cval_nlcc(3, iexp_nlcc)/alpha**4
381 0 : n = coset(2, 2, 0)
382 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
383 0 : n = coset(2, 0, 2)
384 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
385 0 : n = coset(0, 2, 2)
386 0 : pab(n, 1) = 2._dp*cval_nlcc(3, iexp_nlcc)/alpha**4
387 : CASE (4)
388 0 : n = coset(6, 0, 0)
389 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
390 0 : n = coset(0, 6, 0)
391 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
392 0 : n = coset(0, 0, 6)
393 0 : pab(n, 1) = cval_nlcc(4, iexp_nlcc)/alpha**6
394 0 : n = coset(4, 2, 0)
395 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
396 0 : n = coset(4, 0, 2)
397 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
398 0 : n = coset(2, 4, 0)
399 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
400 0 : n = coset(2, 0, 4)
401 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
402 0 : n = coset(0, 4, 2)
403 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
404 0 : n = coset(0, 2, 4)
405 0 : pab(n, 1) = 3._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
406 0 : n = coset(2, 2, 2)
407 0 : pab(n, 1) = 6._dp*cval_nlcc(4, iexp_nlcc)/alpha**6
408 : CASE DEFAULT
409 58 : CPABORT("")
410 : END SELECT
411 : END DO
412 42 : IF (dft_control%nspins == 2) pab = pab*0.5_dp
413 :
414 172 : DO iatom = 1, natom
415 130 : atom_a = atom_list(iatom)
416 130 : ra(:) = pbc(particle_set(atom_a)%r, cell)
417 172 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
418 : ! replicated realspace grid, split the atoms up between procs
419 130 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
420 65 : npme = npme + 1
421 65 : cores(npme) = iatom
422 : END IF
423 : ELSE
424 0 : npme = npme + 1
425 0 : cores(npme) = iatom
426 : END IF
427 : END DO
428 :
429 107 : DO j = 1, npme
430 :
431 65 : iatom = cores(j)
432 65 : atom_a = atom_list(iatom)
433 65 : ra(:) = pbc(particle_set(atom_a)%r, cell)
434 274 : hab(:, 1) = 0.0_dp
435 65 : force_a(:) = 0.0_dp
436 65 : force_b(:) = 0.0_dp
437 65 : IF (use_virial) THEN
438 48 : my_virial_a = 0.0_dp
439 48 : my_virial_b = 0.0_dp
440 : END IF
441 65 : ni = 2*nc - 2
442 :
443 : radius = exp_radius_very_extended(la_min=0, la_max=ni, lb_min=0, lb_max=0, &
444 : ra=ra, rb=ra, rp=ra, &
445 : zetp=1/(2*alpha**2), eps=eps_rho_rspace, &
446 : pab=pab, o1=0, o2=0, & ! without map_consistent
447 65 : prefactor=1.0_dp, cutoff=1.0_dp)
448 :
449 : CALL integrate_pgf_product(ni, 1/(2*alpha**2), 0, &
450 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
451 : rs_v, hab, pab=pab, o1=0, o2=0, &
452 : radius=radius, &
453 : calculate_forces=.TRUE., force_a=force_a, &
454 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
455 65 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
456 :
457 : force(ikind)%gth_nlcc(:, iatom) = &
458 260 : force(ikind)%gth_nlcc(:, iatom) + force_a(:)*rho_rspace%pw_grid%dvol
459 :
460 107 : IF (use_virial) THEN
461 624 : virial%pv_nlcc = virial%pv_nlcc + my_virial_a*rho_rspace%pw_grid%dvol
462 624 : virial%pv_virial = virial%pv_virial + my_virial_a*rho_rspace%pw_grid%dvol
463 : END IF
464 : END DO
465 :
466 86 : DEALLOCATE (hab, pab)
467 :
468 : END DO
469 :
470 : END DO
471 :
472 30 : DEALLOCATE (cores)
473 :
474 30 : CALL timestop(handle)
475 :
476 30 : END SUBROUTINE integrate_rho_nlcc
477 :
478 : ! **************************************************************************************************
479 : !> \brief computes the forces/virial due to the ionic cores with a potential on
480 : !> grid
481 : !> \param v_rspace ...
482 : !> \param qs_env ...
483 : !> \param atecc ...
484 : ! **************************************************************************************************
485 7531 : SUBROUTINE integrate_v_core_rspace(v_rspace, qs_env, atecc)
486 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
487 : TYPE(qs_environment_type), POINTER :: qs_env
488 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atecc
489 :
490 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_core_rspace'
491 :
492 : INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
493 : natom_of_kind, npme
494 7531 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
495 : LOGICAL :: paw_atom, skip_fcore, use_virial
496 : REAL(KIND=dp) :: alpha_core_charge, ccore_charge, &
497 : eps_rho_rspace, radius
498 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
499 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
500 7531 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
501 7531 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
502 : TYPE(atprop_type), POINTER :: atprop
503 : TYPE(cell_type), POINTER :: cell
504 : TYPE(dft_control_type), POINTER :: dft_control
505 7531 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
506 : TYPE(pw_env_type), POINTER :: pw_env
507 7531 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
508 7531 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
509 : TYPE(realspace_grid_type), POINTER :: rs_v
510 : TYPE(virial_type), POINTER :: virial
511 :
512 7531 : CALL timeset(routineN, handle)
513 7531 : NULLIFY (virial, force, atprop, dft_control)
514 :
515 7531 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
516 :
517 : !If gapw, check for gpw kinds
518 7531 : skip_fcore = .FALSE.
519 7531 : IF (dft_control%qs_control%gapw) THEN
520 546 : IF (.NOT. dft_control%qs_control%gapw_control%nopaw_as_gpw) skip_fcore = .TRUE.
521 : END IF
522 :
523 : IF (.NOT. skip_fcore) THEN
524 7043 : NULLIFY (pw_env)
525 7043 : ALLOCATE (cores(1))
526 7043 : ALLOCATE (hab(1, 1))
527 7043 : ALLOCATE (pab(1, 1))
528 :
529 7043 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
530 7043 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
531 :
532 7043 : CALL transfer_pw2rs(rs_v, v_rspace)
533 :
534 : CALL get_qs_env(qs_env=qs_env, &
535 : atomic_kind_set=atomic_kind_set, &
536 : qs_kind_set=qs_kind_set, &
537 : cell=cell, &
538 : dft_control=dft_control, &
539 : particle_set=particle_set, &
540 : pw_env=pw_env, &
541 : force=force, &
542 : virial=virial, &
543 7043 : atprop=atprop)
544 :
545 : ! atomic energy contributions
546 7043 : natom = SIZE(particle_set)
547 7043 : IF (ASSOCIATED(atprop)) THEN
548 7043 : CALL atprop_array_init(atprop%ateb, natom)
549 : END IF
550 :
551 7043 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
552 :
553 7043 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
554 :
555 19413 : DO ikind = 1, SIZE(atomic_kind_set)
556 :
557 12370 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
558 12370 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
559 :
560 12370 : IF (dft_control%qs_control%gapw .AND. paw_atom) CYCLE
561 :
562 : CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha_core_charge, &
563 12306 : ccore_charge=ccore_charge)
564 :
565 12306 : pab(1, 1) = -ccore_charge
566 12306 : IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
567 :
568 12270 : CALL reallocate(cores, 1, natom_of_kind)
569 12270 : npme = 0
570 37111 : cores = 0
571 :
572 37111 : DO iatom = 1, natom_of_kind
573 24841 : atom_a = atom_list(iatom)
574 24841 : ra(:) = pbc(particle_set(atom_a)%r, cell)
575 37111 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
576 : ! replicated realspace grid, split the atoms up between procs
577 24092 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
578 12046 : npme = npme + 1
579 12046 : cores(npme) = iatom
580 : END IF
581 : ELSE
582 749 : npme = npme + 1
583 749 : cores(npme) = iatom
584 : END IF
585 : END DO
586 :
587 56784 : DO j = 1, npme
588 :
589 12795 : iatom = cores(j)
590 12795 : atom_a = atom_list(iatom)
591 12795 : ra(:) = pbc(particle_set(atom_a)%r, cell)
592 12795 : hab(1, 1) = 0.0_dp
593 12795 : force_a(:) = 0.0_dp
594 12795 : force_b(:) = 0.0_dp
595 12795 : IF (use_virial) THEN
596 1600 : my_virial_a = 0.0_dp
597 1600 : my_virial_b = 0.0_dp
598 : END IF
599 :
600 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
601 : ra=ra, rb=ra, rp=ra, &
602 : zetp=alpha_core_charge, eps=eps_rho_rspace, &
603 : pab=pab, o1=0, o2=0, & ! without map_consistent
604 12795 : prefactor=1.0_dp, cutoff=1.0_dp)
605 :
606 : CALL integrate_pgf_product(0, alpha_core_charge, 0, &
607 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
608 : rs_v, hab, pab=pab, o1=0, o2=0, &
609 : radius=radius, &
610 : calculate_forces=.TRUE., force_a=force_a, &
611 : force_b=force_b, use_virial=use_virial, my_virial_a=my_virial_a, &
612 12795 : my_virial_b=my_virial_b, use_subpatch=.TRUE., subpatch_pattern=0)
613 :
614 12795 : IF (ASSOCIATED(force)) THEN
615 50792 : force(ikind)%rho_core(:, iatom) = force(ikind)%rho_core(:, iatom) + force_a(:)
616 : END IF
617 12795 : IF (use_virial) THEN
618 20800 : virial%pv_ehartree = virial%pv_ehartree + my_virial_a
619 20800 : virial%pv_virial = virial%pv_virial + my_virial_a
620 : END IF
621 12795 : IF (ASSOCIATED(atprop)) THEN
622 12795 : atprop%ateb(atom_a) = atprop%ateb(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
623 : END IF
624 25165 : IF (PRESENT(atecc)) THEN
625 47 : atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
626 : END IF
627 :
628 : END DO
629 :
630 : END DO
631 :
632 7043 : DEALLOCATE (hab, pab, cores)
633 :
634 : END IF
635 :
636 7531 : CALL timestop(handle)
637 :
638 7531 : END SUBROUTINE integrate_v_core_rspace
639 :
640 : ! **************************************************************************************************
641 : !> \brief computes the overlap of a set of Gaussians with a potential on grid
642 : !> \param v_rspace ...
643 : !> \param qs_env ...
644 : !> \param alpha ...
645 : !> \param ccore ...
646 : !> \param atecc ...
647 : ! **************************************************************************************************
648 2 : SUBROUTINE integrate_v_gaussian_rspace(v_rspace, qs_env, alpha, ccore, atecc)
649 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
650 : TYPE(qs_environment_type), POINTER :: qs_env
651 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: alpha, ccore
652 : REAL(KIND=dp), DIMENSION(:) :: atecc
653 :
654 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_gaussian_rspace'
655 :
656 : INTEGER :: atom_a, handle, iatom, ikind, j, natom, &
657 : natom_of_kind, npme
658 2 : INTEGER, DIMENSION(:), POINTER :: atom_list, cores
659 : REAL(KIND=dp) :: alpha_core_charge, eps_rho_rspace, radius
660 : REAL(KIND=dp), DIMENSION(3) :: ra
661 2 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab
662 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
663 : TYPE(cell_type), POINTER :: cell
664 : TYPE(dft_control_type), POINTER :: dft_control
665 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
666 : TYPE(pw_env_type), POINTER :: pw_env
667 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
668 : TYPE(realspace_grid_type), POINTER :: rs_v
669 :
670 2 : CALL timeset(routineN, handle)
671 :
672 2 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
673 :
674 : !If gapw, check for gpw kinds
675 2 : CPASSERT(.NOT. dft_control%qs_control%gapw)
676 :
677 2 : NULLIFY (pw_env)
678 2 : ALLOCATE (cores(1))
679 2 : ALLOCATE (hab(1, 1))
680 2 : ALLOCATE (pab(1, 1))
681 :
682 2 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
683 2 : CALL pw_env_get(pw_env=pw_env, auxbas_rs_grid=rs_v)
684 :
685 2 : CALL transfer_pw2rs(rs_v, v_rspace)
686 :
687 : CALL get_qs_env(qs_env=qs_env, &
688 : atomic_kind_set=atomic_kind_set, &
689 : qs_kind_set=qs_kind_set, &
690 : cell=cell, &
691 : dft_control=dft_control, &
692 : particle_set=particle_set, &
693 2 : pw_env=pw_env)
694 :
695 : ! atomic energy contributions
696 2 : natom = SIZE(particle_set)
697 2 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
698 :
699 6 : DO ikind = 1, SIZE(atomic_kind_set)
700 :
701 4 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
702 4 : pab(1, 1) = -ccore(ikind)
703 4 : alpha_core_charge = alpha(ikind)
704 4 : IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE
705 :
706 4 : CALL reallocate(cores, 1, natom_of_kind)
707 4 : npme = 0
708 10 : cores = 0
709 :
710 10 : DO iatom = 1, natom_of_kind
711 6 : atom_a = atom_list(iatom)
712 6 : ra(:) = pbc(particle_set(atom_a)%r, cell)
713 10 : IF (rs_v%desc%parallel .AND. .NOT. rs_v%desc%distributed) THEN
714 : ! replicated realspace grid, split the atoms up between procs
715 6 : IF (MODULO(iatom, rs_v%desc%group_size) == rs_v%desc%my_pos) THEN
716 3 : npme = npme + 1
717 3 : cores(npme) = iatom
718 : END IF
719 : ELSE
720 0 : npme = npme + 1
721 0 : cores(npme) = iatom
722 : END IF
723 : END DO
724 :
725 13 : DO j = 1, npme
726 :
727 3 : iatom = cores(j)
728 3 : atom_a = atom_list(iatom)
729 3 : ra(:) = pbc(particle_set(atom_a)%r, cell)
730 3 : hab(1, 1) = 0.0_dp
731 :
732 : radius = exp_radius_very_extended(la_min=0, la_max=0, lb_min=0, lb_max=0, &
733 : ra=ra, rb=ra, rp=ra, &
734 : zetp=alpha_core_charge, eps=eps_rho_rspace, &
735 : pab=pab, o1=0, o2=0, & ! without map_consistent
736 3 : prefactor=1.0_dp, cutoff=1.0_dp)
737 :
738 : CALL integrate_pgf_product(0, alpha_core_charge, 0, &
739 : 0, 0.0_dp, 0, ra, [0.0_dp, 0.0_dp, 0.0_dp], &
740 : rs_v, hab, pab=pab, o1=0, o2=0, &
741 : radius=radius, calculate_forces=.FALSE., &
742 3 : use_subpatch=.TRUE., subpatch_pattern=0)
743 7 : atecc(atom_a) = atecc(atom_a) + 0.5_dp*hab(1, 1)*pab(1, 1)
744 :
745 : END DO
746 :
747 : END DO
748 :
749 2 : DEALLOCATE (hab, pab, cores)
750 :
751 2 : CALL timestop(handle)
752 :
753 2 : END SUBROUTINE integrate_v_gaussian_rspace
754 : ! **************************************************************************************************
755 : !> \brief computes integrals of product of v_rspace times a one-center function
756 : !> required for LRIGPW
757 : !> \param v_rspace ...
758 : !> \param qs_env ...
759 : !> \param int_res ...
760 : !> \param calculate_forces ...
761 : !> \param basis_type ...
762 : !> \param atomlist ...
763 : !> \author Dorothea Golze
764 : ! **************************************************************************************************
765 822 : SUBROUTINE integrate_v_rspace_one_center(v_rspace, qs_env, int_res, &
766 822 : calculate_forces, basis_type, atomlist)
767 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
768 : TYPE(qs_environment_type), POINTER :: qs_env
769 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: int_res
770 : LOGICAL, INTENT(IN) :: calculate_forces
771 : CHARACTER(len=*), INTENT(IN) :: basis_type
772 : INTEGER, DIMENSION(:), OPTIONAL :: atomlist
773 :
774 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_one_center'
775 :
776 : INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, m1, maxco, &
777 : maxsgf_set, my_pos, na1, natom_of_kind, ncoa, nkind, nseta, offset, sgfa
778 822 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
779 822 : nsgf_seta
780 822 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
781 : LOGICAL :: use_virial
782 822 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: map_it
783 : REAL(KIND=dp) :: eps_rho_rspace, radius
784 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
785 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
786 822 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
787 822 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab, rpgfa, sphi_a, work_f, work_i, &
788 822 : zeta
789 822 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
790 : TYPE(cell_type), POINTER :: cell
791 : TYPE(dft_control_type), POINTER :: dft_control
792 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
793 : TYPE(gto_basis_set_type), POINTER :: lri_basis_set
794 822 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
795 : TYPE(pw_env_type), POINTER :: pw_env
796 822 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
797 822 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
798 : TYPE(realspace_grid_type), POINTER :: rs_grid
799 : TYPE(virial_type), POINTER :: virial
800 :
801 822 : CALL timeset(routineN, handle)
802 :
803 822 : NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
804 822 : first_sgfa, gridlevel_info, hab, la_max, la_min, lri_basis_set, &
805 822 : npgfa, nsgf_seta, pab, particle_set, pw_env, rpgfa, &
806 822 : rs_grid, rs_v, virial, set_radius_a, sphi_a, work_f, &
807 822 : work_i, zeta)
808 :
809 822 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
810 :
811 822 : CALL pw_env_get(pw_env, rs_grids=rs_v)
812 4074 : DO i = 1, SIZE(rs_v)
813 4074 : CALL rs_grid_zero(rs_v(i))
814 : END DO
815 :
816 822 : gridlevel_info => pw_env%gridlevel_info
817 :
818 822 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
819 :
820 : CALL get_qs_env(qs_env=qs_env, &
821 : atomic_kind_set=atomic_kind_set, &
822 : qs_kind_set=qs_kind_set, &
823 : cell=cell, &
824 : dft_control=dft_control, &
825 : nkind=nkind, &
826 : particle_set=particle_set, &
827 : pw_env=pw_env, &
828 822 : virial=virial)
829 :
830 822 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
831 :
832 822 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
833 :
834 822 : offset = 0
835 822 : my_pos = v_rspace%pw_grid%para%group%mepos
836 822 : group_size = v_rspace%pw_grid%para%group%num_pe
837 :
838 2448 : DO ikind = 1, nkind
839 :
840 1626 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
841 1626 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=lri_basis_set, basis_type=basis_type)
842 : CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
843 : first_sgf=first_sgfa, &
844 : lmax=la_max, &
845 : lmin=la_min, &
846 : maxco=maxco, &
847 : maxsgf_set=maxsgf_set, &
848 : npgf=npgfa, &
849 : nset=nseta, &
850 : nsgf_set=nsgf_seta, &
851 : pgf_radius=rpgfa, &
852 : set_radius=set_radius_a, &
853 : sphi=sphi_a, &
854 1626 : zet=zeta)
855 :
856 6504 : ALLOCATE (hab(maxco, 1), pab(maxco, 1))
857 43850 : hab = 0._dp
858 42224 : pab(:, 1) = 0._dp
859 :
860 4746 : DO iatom = 1, natom_of_kind
861 :
862 3120 : atom_a = atom_list(iatom)
863 3120 : IF (PRESENT(atomlist)) THEN
864 400 : IF (atomlist(atom_a) == 0) CYCLE
865 : END IF
866 2920 : ra(:) = pbc(particle_set(atom_a)%r, cell)
867 2920 : force_a(:) = 0._dp
868 2920 : force_b(:) = 0._dp
869 2920 : my_virial_a(:, :) = 0._dp
870 2920 : my_virial_b(:, :) = 0._dp
871 :
872 42868 : m1 = MAXVAL(npgfa(1:nseta))
873 8760 : ALLOCATE (map_it(m1))
874 :
875 42868 : DO iset = 1, nseta
876 : !
877 84320 : map_it = .FALSE.
878 84128 : DO ipgf = 1, npgfa(iset)
879 44180 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
880 44180 : rs_grid => rs_v(igrid_level)
881 84128 : map_it(ipgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
882 : END DO
883 39948 : offset = offset + 1
884 : !
885 64958 : IF (ANY(map_it(1:npgfa(iset)))) THEN
886 19974 : sgfa = first_sgfa(1, iset)
887 19974 : ncoa = npgfa(iset)*ncoset(la_max(iset))
888 424846 : hab(:, 1) = 0._dp
889 59922 : ALLOCATE (work_i(nsgf_seta(iset), 1))
890 192924 : work_i = 0.0_dp
891 :
892 : ! get fit coefficients for forces
893 19974 : IF (calculate_forces) THEN
894 895 : m1 = sgfa + nsgf_seta(iset) - 1
895 1790 : ALLOCATE (work_f(nsgf_seta(iset), 1))
896 6465 : work_f(1:nsgf_seta(iset), 1) = int_res(ikind)%acoef(iatom, sgfa:m1)
897 : CALL dgemm("N", "N", ncoa, 1, nsgf_seta(iset), 1.0_dp, sphi_a(1, sgfa), &
898 : SIZE(sphi_a, 1), work_f(1, 1), SIZE(work_f, 1), 0.0_dp, pab(1, 1), &
899 895 : SIZE(pab, 1))
900 895 : DEALLOCATE (work_f)
901 : END IF
902 :
903 42064 : DO ipgf = 1, npgfa(iset)
904 22090 : na1 = (ipgf - 1)*ncoset(la_max(iset))
905 22090 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
906 22090 : rs_grid => rs_v(igrid_level)
907 :
908 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
909 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
910 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
911 22090 : prefactor=1.0_dp, cutoff=1.0_dp)
912 :
913 42064 : IF (map_it(ipgf)) THEN
914 22090 : IF (.NOT. calculate_forces) THEN
915 : CALL integrate_pgf_product(la_max=la_max(iset), &
916 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
917 : lb_max=0, zetb=0.0_dp, lb_min=0, &
918 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
919 : rsgrid=rs_grid, &
920 : hab=hab, o1=na1, o2=0, radius=radius, &
921 20831 : calculate_forces=calculate_forces)
922 : ELSE
923 : CALL integrate_pgf_product(la_max=la_max(iset), &
924 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
925 : lb_max=0, zetb=0.0_dp, lb_min=0, &
926 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
927 : rsgrid=rs_grid, &
928 : hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
929 : calculate_forces=calculate_forces, &
930 : force_a=force_a, force_b=force_b, &
931 : use_virial=use_virial, &
932 1259 : my_virial_a=my_virial_a, my_virial_b=my_virial_b)
933 : END IF
934 : END IF
935 : END DO
936 : ! contract hab
937 : CALL dgemm("T", "N", nsgf_seta(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
938 19974 : SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work_i(1, 1), SIZE(work_i, 1))
939 :
940 : int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) = &
941 325926 : int_res(ikind)%v_int(iatom, sgfa:sgfa - 1 + nsgf_seta(iset)) + work_i(1:nsgf_seta(iset), 1)
942 19974 : DEALLOCATE (work_i)
943 : END IF
944 : END DO
945 : !
946 2920 : IF (calculate_forces) THEN
947 568 : int_res(ikind)%v_dfdr(iatom, :) = int_res(ikind)%v_dfdr(iatom, :) + force_a(:)
948 142 : IF (use_virial) THEN
949 156 : virial%pv_lrigpw = virial%pv_lrigpw + my_virial_a
950 156 : virial%pv_virial = virial%pv_virial + my_virial_a
951 : END IF
952 : END IF
953 :
954 4746 : DEALLOCATE (map_it)
955 :
956 : END DO
957 :
958 5700 : DEALLOCATE (hab, pab)
959 : END DO
960 :
961 822 : CALL timestop(handle)
962 :
963 1644 : END SUBROUTINE integrate_v_rspace_one_center
964 :
965 : ! **************************************************************************************************
966 : !> \brief computes integrals of product of v_rspace times the diagonal block basis functions
967 : !> required for LRIGPW with exact 1c terms
968 : !> \param v_rspace ...
969 : !> \param ksmat ...
970 : !> \param pmat ...
971 : !> \param qs_env ...
972 : !> \param calculate_forces ...
973 : !> \param basis_type ...
974 : !> \author JGH
975 : ! **************************************************************************************************
976 36 : SUBROUTINE integrate_v_rspace_diagonal(v_rspace, ksmat, pmat, qs_env, calculate_forces, basis_type)
977 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
978 : TYPE(dbcsr_type), INTENT(INOUT) :: ksmat, pmat
979 : TYPE(qs_environment_type), POINTER :: qs_env
980 : LOGICAL, INTENT(IN) :: calculate_forces
981 : CHARACTER(len=*), INTENT(IN) :: basis_type
982 :
983 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_v_rspace_diagonal'
984 :
985 : INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
986 : m1, maxco, maxsgf_set, my_pos, na1, na2, natom_of_kind, nb1, nb2, ncoa, ncob, nkind, &
987 : nseta, nsgfa, offset, sgfa, sgfb
988 36 : INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, &
989 36 : nsgf_seta
990 36 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
991 : LOGICAL :: found, use_virial
992 36 : LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
993 : REAL(KIND=dp) :: eps_rho_rspace, radius, zetp
994 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
995 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
996 36 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
997 36 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, hab, hmat, p_block, pab, pblk, &
998 36 : rpgfa, sphi_a, work, zeta
999 36 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1000 : TYPE(cell_type), POINTER :: cell
1001 : TYPE(dft_control_type), POINTER :: dft_control
1002 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1003 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1004 : TYPE(mp_para_env_type), POINTER :: para_env
1005 36 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1006 : TYPE(pw_env_type), POINTER :: pw_env
1007 36 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1008 36 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1009 36 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1010 : TYPE(realspace_grid_type), POINTER :: rs_grid
1011 : TYPE(virial_type), POINTER :: virial
1012 :
1013 36 : CALL timeset(routineN, handle)
1014 :
1015 36 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1016 36 : CALL pw_env_get(pw_env, rs_grids=rs_v)
1017 36 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1018 :
1019 36 : gridlevel_info => pw_env%gridlevel_info
1020 :
1021 : CALL get_qs_env(qs_env=qs_env, &
1022 : atomic_kind_set=atomic_kind_set, &
1023 : qs_kind_set=qs_kind_set, &
1024 : cell=cell, &
1025 : dft_control=dft_control, &
1026 : nkind=nkind, &
1027 : particle_set=particle_set, &
1028 : force=force, &
1029 : virial=virial, &
1030 36 : para_env=para_env)
1031 :
1032 36 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1033 36 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1034 :
1035 36 : offset = 0
1036 36 : my_pos = v_rspace%pw_grid%para%group%mepos
1037 36 : group_size = v_rspace%pw_grid%para%group%num_pe
1038 :
1039 108 : DO ikind = 1, nkind
1040 :
1041 72 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
1042 72 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1043 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1044 : lmax=la_max, lmin=la_min, maxco=maxco, maxsgf_set=maxsgf_set, &
1045 : npgf=npgfa, nset=nseta, nsgf_set=nsgf_seta, nsgf=nsgfa, &
1046 : first_sgf=first_sgfa, pgf_radius=rpgfa, set_radius=set_radius_a, &
1047 72 : sphi=sphi_a, zet=zeta)
1048 :
1049 720 : ALLOCATE (hab(maxco, maxco), work(maxco, maxsgf_set), hmat(nsgfa, nsgfa))
1050 72 : IF (calculate_forces) ALLOCATE (pab(maxco, maxco), pblk(nsgfa, nsgfa))
1051 :
1052 288 : DO iatom = 1, natom_of_kind
1053 216 : atom_a = atom_list(iatom)
1054 216 : ra(:) = pbc(particle_set(atom_a)%r, cell)
1055 17640 : hmat = 0.0_dp
1056 216 : IF (calculate_forces) THEN
1057 0 : CALL dbcsr_get_block_p(matrix=pmat, row=atom_a, col=atom_a, BLOCK=p_block, found=found)
1058 0 : IF (found) THEN
1059 0 : pblk(1:nsgfa, 1:nsgfa) = p_block(1:nsgfa, 1:nsgfa)
1060 : ELSE
1061 0 : pblk = 0.0_dp
1062 : END IF
1063 0 : CALL para_env%sum(pblk)
1064 0 : force_a(:) = 0._dp
1065 0 : force_b(:) = 0._dp
1066 0 : IF (use_virial) THEN
1067 0 : my_virial_a = 0.0_dp
1068 0 : my_virial_b = 0.0_dp
1069 : END IF
1070 : END IF
1071 432 : m1 = MAXVAL(npgfa(1:nseta))
1072 864 : ALLOCATE (map_it2(m1, m1))
1073 432 : DO iset = 1, nseta
1074 216 : sgfa = first_sgfa(1, iset)
1075 216 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1076 648 : DO jset = 1, nseta
1077 216 : sgfb = first_sgfa(1, jset)
1078 216 : ncob = npgfa(jset)*ncoset(la_max(jset))
1079 : !
1080 12312 : map_it2 = .FALSE.
1081 1728 : DO ipgf = 1, npgfa(iset)
1082 12312 : DO jpgf = 1, npgfa(jset)
1083 10584 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1084 10584 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1085 10584 : rs_grid => rs_v(igrid_level)
1086 12096 : map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
1087 : END DO
1088 : END DO
1089 216 : offset = offset + 1
1090 : !
1091 6480 : IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
1092 237492 : hab(:, :) = 0._dp
1093 108 : IF (calculate_forces) THEN
1094 : CALL dgemm("N", "N", ncoa, nsgf_seta(jset), nsgf_seta(iset), &
1095 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1096 : pblk(sgfa, sgfb), SIZE(pblk, 1), &
1097 0 : 0.0_dp, work(1, 1), SIZE(work, 1))
1098 : CALL dgemm("N", "T", ncoa, ncob, nsgf_seta(jset), &
1099 : 1.0_dp, work(1, 1), SIZE(work, 1), &
1100 : sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1101 0 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
1102 : END IF
1103 :
1104 864 : DO ipgf = 1, npgfa(iset)
1105 756 : na1 = (ipgf - 1)*ncoset(la_max(iset))
1106 756 : na2 = ipgf*ncoset(la_max(iset))
1107 6156 : DO jpgf = 1, npgfa(jset)
1108 5292 : nb1 = (jpgf - 1)*ncoset(la_max(jset))
1109 5292 : nb2 = jpgf*ncoset(la_max(jset))
1110 5292 : zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
1111 5292 : igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
1112 5292 : rs_grid => rs_v(igrid_level)
1113 :
1114 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1115 : lb_min=la_min(jset), lb_max=la_max(jset), &
1116 : ra=ra, rb=ra, rp=ra, &
1117 : zetp=zetp, eps=eps_rho_rspace, &
1118 5292 : prefactor=1.0_dp, cutoff=1.0_dp)
1119 :
1120 6048 : IF (map_it2(ipgf, jpgf)) THEN
1121 5292 : IF (calculate_forces) THEN
1122 : CALL integrate_pgf_product( &
1123 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
1124 : la_max(jset), zeta(jpgf, jset), la_min(jset), &
1125 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
1126 : rsgrid=rs_v(igrid_level), &
1127 : hab=hab, pab=pab, o1=na1, o2=nb1, &
1128 : radius=radius, &
1129 : calculate_forces=.TRUE., &
1130 : force_a=force_a, force_b=force_b, &
1131 0 : use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1132 : ELSE
1133 : CALL integrate_pgf_product( &
1134 : la_max(iset), zeta(ipgf, iset), la_min(iset), &
1135 : la_max(jset), zeta(jpgf, jset), la_min(jset), &
1136 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
1137 : rsgrid=rs_v(igrid_level), &
1138 : hab=hab, o1=na1, o2=nb1, &
1139 : radius=radius, &
1140 5292 : calculate_forces=.FALSE.)
1141 : END IF
1142 : END IF
1143 : END DO
1144 : END DO
1145 : ! contract hab
1146 : CALL dgemm("N", "N", ncoa, nsgf_seta(jset), ncob, &
1147 : 1.0_dp, hab(1, 1), SIZE(hab, 1), &
1148 : sphi_a(1, sgfb), SIZE(sphi_a, 1), &
1149 108 : 0.0_dp, work(1, 1), SIZE(work, 1))
1150 : CALL dgemm("T", "N", nsgf_seta(iset), nsgf_seta(jset), ncoa, &
1151 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1152 : work(1, 1), SIZE(work, 1), &
1153 108 : 1.0_dp, hmat(sgfa, sgfb), SIZE(hmat, 1))
1154 : END IF
1155 : END DO
1156 : END DO
1157 216 : DEALLOCATE (map_it2)
1158 : ! update KS matrix
1159 35064 : CALL para_env%sum(hmat)
1160 216 : CALL dbcsr_get_block_p(matrix=ksmat, row=atom_a, col=atom_a, BLOCK=h_block, found=found)
1161 216 : IF (found) THEN
1162 17532 : h_block(1:nsgfa, 1:nsgfa) = h_block(1:nsgfa, 1:nsgfa) + hmat(1:nsgfa, 1:nsgfa)
1163 : END IF
1164 504 : IF (calculate_forces) THEN
1165 0 : force(ikind)%rho_elec(:, iatom) = force(ikind)%rho_elec(:, iatom) + 2.0_dp*force_a(:)
1166 0 : IF (use_virial) THEN
1167 : IF (use_virial .AND. calculate_forces) THEN
1168 0 : virial%pv_lrigpw = virial%pv_lrigpw + 2.0_dp*my_virial_a
1169 0 : virial%pv_virial = virial%pv_virial + 2.0_dp*my_virial_a
1170 : END IF
1171 : END IF
1172 : END IF
1173 : END DO
1174 72 : DEALLOCATE (hab, work, hmat)
1175 252 : IF (calculate_forces) DEALLOCATE (pab, pblk)
1176 : END DO
1177 :
1178 36 : CALL timestop(handle)
1179 :
1180 72 : END SUBROUTINE integrate_v_rspace_diagonal
1181 :
1182 : ! **************************************************************************************************
1183 : !> \brief computes integrals of product of v_rspace times a basis function (vector function)
1184 : !> and possible forces
1185 : !> \param qs_env ...
1186 : !> \param v_rspace ...
1187 : !> \param f_coef ...
1188 : !> \param f_integral ...
1189 : !> \param calculate_forces ...
1190 : !> \param basis_type ...
1191 : !> \author JGH [8.2024]
1192 : ! **************************************************************************************************
1193 4 : SUBROUTINE integrate_function(qs_env, v_rspace, f_coef, f_integral, calculate_forces, basis_type)
1194 : TYPE(qs_environment_type), POINTER :: qs_env
1195 : TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
1196 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: f_coef
1197 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: f_integral
1198 : LOGICAL, INTENT(IN) :: calculate_forces
1199 : CHARACTER(len=*), INTENT(IN) :: basis_type
1200 :
1201 : CHARACTER(len=*), PARAMETER :: routineN = 'integrate_function'
1202 :
1203 : INTEGER :: atom_a, group_size, handle, i, iatom, igrid_level, ikind, ipgf, iset, maxco, &
1204 : maxsgf_set, my_pos, na1, natom, ncoa, nkind, nseta, offset, sgfa
1205 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
1206 4 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
1207 4 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
1208 : LOGICAL :: use_virial
1209 : REAL(KIND=dp) :: eps_rho_rspace, radius
1210 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
1211 : REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
1212 4 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, pab, sphi_a, work, zeta
1213 4 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1214 : TYPE(cell_type), POINTER :: cell
1215 : TYPE(dft_control_type), POINTER :: dft_control
1216 : TYPE(gridlevel_info_type), POINTER :: gridlevel_info
1217 : TYPE(gto_basis_set_type), POINTER :: orb_basis_set
1218 4 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1219 : TYPE(pw_env_type), POINTER :: pw_env
1220 4 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1221 4 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1222 4 : TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_v
1223 : TYPE(realspace_grid_type), POINTER :: rs_grid
1224 : TYPE(virial_type), POINTER :: virial
1225 :
1226 4 : CALL timeset(routineN, handle)
1227 :
1228 4 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env)
1229 4 : gridlevel_info => pw_env%gridlevel_info
1230 :
1231 4 : CALL pw_env_get(pw_env, rs_grids=rs_v)
1232 4 : CALL potential_pw2rs(rs_v, v_rspace, pw_env)
1233 :
1234 : CALL get_qs_env(qs_env=qs_env, &
1235 : atomic_kind_set=atomic_kind_set, &
1236 : qs_kind_set=qs_kind_set, &
1237 : force=force, &
1238 : cell=cell, &
1239 : dft_control=dft_control, &
1240 : nkind=nkind, &
1241 : natom=natom, &
1242 : particle_set=particle_set, &
1243 : pw_env=pw_env, &
1244 4 : virial=virial)
1245 4 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1246 :
1247 4 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
1248 4 : IF (use_virial) THEN
1249 0 : CPABORT("Virial NYA")
1250 : END IF
1251 :
1252 4 : eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
1253 :
1254 : CALL get_qs_kind_set(qs_kind_set, &
1255 4 : maxco=maxco, maxsgf_set=maxsgf_set, basis_type=basis_type)
1256 20 : ALLOCATE (hab(maxco, 1), pab(maxco, 1), work(maxco, 1))
1257 :
1258 4 : offset = 0
1259 4 : my_pos = v_rspace%pw_grid%para%group%mepos
1260 4 : group_size = v_rspace%pw_grid%para%group%num_pe
1261 :
1262 20 : DO iatom = 1, natom
1263 16 : ikind = particle_set(iatom)%atomic_kind%kind_number
1264 16 : atom_a = atom_of_kind(iatom)
1265 16 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type=basis_type)
1266 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
1267 : first_sgf=first_sgfa, &
1268 : lmax=la_max, &
1269 : lmin=la_min, &
1270 : npgf=npgfa, &
1271 : nset=nseta, &
1272 : nsgf_set=nsgfa, &
1273 : sphi=sphi_a, &
1274 16 : zet=zeta)
1275 16 : ra(:) = pbc(particle_set(iatom)%r, cell)
1276 :
1277 16 : force_a(:) = 0._dp
1278 16 : force_b(:) = 0._dp
1279 16 : my_virial_a(:, :) = 0._dp
1280 16 : my_virial_b(:, :) = 0._dp
1281 :
1282 32 : DO iset = 1, nseta
1283 :
1284 16 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1285 16 : sgfa = first_sgfa(1, iset)
1286 :
1287 144 : hab = 0._dp
1288 144 : pab = 0._dp
1289 :
1290 120 : DO i = 1, nsgfa(iset)
1291 120 : work(i, 1) = f_coef(offset + i)
1292 : END DO
1293 :
1294 : CALL dgemm("N", "N", ncoa, 1, nsgfa(iset), &
1295 : 1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
1296 : work(1, 1), SIZE(work, 1), &
1297 16 : 0.0_dp, pab(1, 1), SIZE(pab, 1))
1298 :
1299 120 : DO ipgf = 1, npgfa(iset)
1300 :
1301 104 : na1 = (ipgf - 1)*ncoset(la_max(iset))
1302 :
1303 104 : igrid_level = gaussian_gridlevel(gridlevel_info, zeta(ipgf, iset))
1304 104 : rs_grid => rs_v(igrid_level)
1305 :
1306 120 : IF (map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)) THEN
1307 : radius = exp_radius_very_extended(la_min=la_min(iset), la_max=la_max(iset), &
1308 : lb_min=0, lb_max=0, ra=ra, rb=ra, rp=ra, &
1309 : zetp=zeta(ipgf, iset), eps=eps_rho_rspace, &
1310 52 : prefactor=1.0_dp, cutoff=1.0_dp)
1311 :
1312 52 : IF (calculate_forces) THEN
1313 : CALL integrate_pgf_product(la_max=la_max(iset), &
1314 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
1315 : lb_max=0, zetb=0.0_dp, lb_min=0, &
1316 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
1317 : rsgrid=rs_grid, &
1318 : hab=hab, pab=pab, o1=na1, o2=0, radius=radius, &
1319 : calculate_forces=.TRUE., &
1320 : force_a=force_a, force_b=force_b, &
1321 : use_virial=use_virial, &
1322 52 : my_virial_a=my_virial_a, my_virial_b=my_virial_b)
1323 : ELSE
1324 : CALL integrate_pgf_product(la_max=la_max(iset), &
1325 : zeta=zeta(ipgf, iset), la_min=la_min(iset), &
1326 : lb_max=0, zetb=0.0_dp, lb_min=0, &
1327 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], &
1328 : rsgrid=rs_grid, &
1329 : hab=hab, o1=na1, o2=0, radius=radius, &
1330 0 : calculate_forces=.FALSE.)
1331 : END IF
1332 :
1333 : END IF
1334 :
1335 : END DO
1336 : !
1337 : CALL dgemm("T", "N", nsgfa(iset), 1, ncoa, 1.0_dp, sphi_a(1, sgfa), &
1338 16 : SIZE(sphi_a, 1), hab(1, 1), SIZE(hab, 1), 0.0_dp, work(1, 1), SIZE(work, 1))
1339 120 : DO i = 1, nsgfa(iset)
1340 120 : f_integral(offset + i) = work(i, 1)
1341 : END DO
1342 :
1343 32 : offset = offset + nsgfa(iset)
1344 :
1345 : END DO
1346 :
1347 36 : IF (calculate_forces) THEN
1348 64 : force(ikind)%rho_elec(:, atom_a) = force(ikind)%rho_elec(:, atom_a) + force_a(:)
1349 16 : IF (use_virial) THEN
1350 0 : virial%pv_virial = virial%pv_virial + my_virial_a
1351 : END IF
1352 : END IF
1353 :
1354 : END DO
1355 :
1356 4 : DEALLOCATE (hab, pab, work)
1357 :
1358 4 : CALL timestop(handle)
1359 :
1360 12 : END SUBROUTINE integrate_function
1361 :
1362 : END MODULE qs_integrate_potential_single
|
